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C/^ . Abstract 



We analyse metapopulation dynamics in terms of an individual-based, stochastic 
model of a finite metapopulation. We suggest a new approach, using the number 
of patches in the population as a large parameter. This approach does not require 
that the number of individuals per patch is large, neither is it necessary to assume a 
time-scale separation between local population dynamics and migration. Our ap- 
proach makes it possible to accurately describe the dynamics of metapopulations 
consisting of many small patches. We focus on metapopulations on the brink of 
extinction. We estimate the time to extinction and describe the most likely path 
to extinction. We find that the logarithm of the time to extinction is proportional 
to the product of two vectors, a vector characterising the distribution of patch 
population sizes in the quasi-steady state, and a vector - related to Fisher's repro- 
duction vector - that quantifies the sensitivity of the quasi-steady state distribu- 
tion to demographic fluctuations. We compare our analytical results to stochastic 
simulations of the model, and discuss the range of validity of the analytical ex- 
pressions. By identifying fast and slow degrees of freedom in the metapopulation 
dynamics, we show that the dynamics of large metapopulations close to extinc- 
tion is approximately described by a deterministic equation originally proposed 
by Levins (1969). We were able to compute the rates in Levins' equation in terms 
of the parameters of our stochastic, individual-based model. It turns out, however, 
that the interpretation of the dynamical variable depends strongly on the intrinsic 
growth rate and carrying capacity of the patches. Only when the growth rate and 
the carrying capacity are large does the slow variable correspond to the number of 
patches, as envisaged by Levins. Last but not least, we discuss how our findings 
relate to other, widely used metapopulation models. 
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1. Introduction 



The habitats of animal populations are often geographically divided into many 
small patches, either because of human interference or because natural habitats are 
patchy. Understanding the dynamics of such populations is a problem of great the- 
oretical and practical interest. Isolated small patches are often extinction prone, 
for example because of inbreed i ng in combination with demographic and envi- 
ronmental stochasticity (|Hanskil,ll999|). When the patches are connected into a 
network by migration, local populations may still be prone to extinction, but the 
whole population may persist because empty patches are re-colonised by migrants 
from surrounding occupied patches. Moreover, when the migration rate is suffi- 
ciently large, local populations can be stabilised by an inflow of immigrants (the 
rescue eff"ect). Given a model of such a population, the main concern is usually 
to find out how the different parameter values in the model affect the growth and 
persistence of the metapopulation. 

In Levins' model, the metapop ulation dynamics is simplified by treating each 
patch as either occupied or empty (lLevinsl ll969). The rates of patch colonisation 
and extinction are expressed as functions of the total fraction of occupied patches 
in the population. In its simplest form Levins' model describes the time change 
of the fraction Q of occupied patches: 



at 



(1) 



Here c is the rate of successful colonisation of an empty patch, and e is the rate at 
which occupied patches turn empty (the rate of patch extinction). How these rates 
are related to the life-history parameters determining the stochastic, individual- 
based metapopulation dynamics is not explicitly known. In spatially explicit mod- 
els the parameters are chose n to be functions of for example th e number of occu- 
pied neighbouring patches (Hui and Lil . l2004l : IRov et al.l . |2008[) . Eq. ([T]) is gener- 
ally motivated by assuming a separation of time scales between colonisation and 
extinction of patc hes on the one hand, and the population dynamics within patches 



on the o t her ha nd (|LevinsL 1 19691 : iHanski and GyllenbergL 1 19931 : iLande et al.L 11998 
Etiennel |2002[) . Many of the mechanisms that can be observed in natural pop 
ulations (e.g. AUee and rescue effects) can be represented by suitable modi 
fications of Eq. ([T]), and the qualitative behaviour of such models is well un- 



derstood (Hanski and Ovaskainen. 


2000; 


Etienne. 


2000; 


Hardine and McNamaraL 


2002; 


Zhou and Wans . 


2004). See also 


Zhou et al.l ( 


2004): 


Tavlor and Hastings 


(2005), and 


Martcheva and Bolkei 


(2007 


). 
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Models that describe the metapopulation dynamics in terms of the fraction of 
occupied patches have no explicit connection to the population dynamics within 
each patch. In other words, it is not apparent how individual b irths, deaths, a nd 
migration events translate into colonisation and extinction rates (lEtiennel |2002|) . 
An alternative line of analysis focuses on the populatio n dynamics within a 



single patch 0rechsler and Wisseill997l : lLopez and Pfisteil . l200ll : lNewman et al. 



20041). The advantage of this approach is that birth and death processes, emigra- 



tion, and immigration can be modelled explicitly in terms of the number of in- 
dividuals in the patch. This makes the model much more immediate in terms of 
biological processes (such as density dependence and AUee effects), but the rest 
of the metapopulation is reduced to a background source of immigrants (assumed 
to be stationary) and is not explicitly modelled. 

Several authors have attempted to bridge the gap between deta iled loca l popu- 
lation dynamics and the dynamics at the overall population level. Keeling (2002j) 
estimated rates of colonisation and extinction events in Levins' equation from 
individual-based simulations. This method rests on the assumption that Eq. ([T]) is 
an appropriate des cription of the metapopulation dynamics. 

Higgins ( 2009 ^ has investigated the extinction risk of metapopulations subject 



to strong dispersal, that is, in a limit where migration is faster than the local popu- 
lation dynamics within patches. As pointed o ut above, the co nverse is commonly 
assumed in writing Eq. ([T])- Within his model, HigginsH 20091) addressed the ques- 
tion of how the degree of fragmentation of the metapopulation affects its risk of 
ext inction. 

IChessonI (1 1981', Il984l), iHanski and GvUenbergl [l99^ Casagrandi and Gattd 
([l999, 2002), Nachman (2000), Metz and Gyllenberg (20o3) and others have anal- 
ysed the equilibrium states and possible persistence of metapopulations in models 
consisting of infin itely many patches with local dynamics coupled via a migration 
pool (reviewed in iMassol et al.L 120090. The assumption that there are infinitely 
many patches is crucial in these studies, as it makes it possible to pose the ques- 
tion under which circumstances the metapopulation persists ad infinitum, that is, 
for which choice of parameters the metapopulation dynamics reaches a non-trivial 
stable equilibrium. For finite metapopulations, by contrast, there is no stable equi- 
librium corresponding to per sistence (for a review of st ochastic extinctions in bi- 
ological populations, see e.g. lOvaskainen and Meersonll2010 ). As is well known, 
finite populations must eventually become extinct (unless they continue to grow). 
This fact is r eferred to as the 'merciless dichotomy of population dynamics' by 



JagersI (| 1992b . Moreover, it remains unclear how the equations employed in these 



studies are related to Eq. ([T]). 
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In this paper we characterise the dynamics of metapopulations with a large, 
but finite, number of patches. We derive, from first principles, how the stochastic 
metapopulation dynamics determines the distribution of individuals over patches. 
The only critical assumption in our derivation is that the number of patches is 
large enough. Especially, it is not necessary to assume that the typical number of 
individuals per patch is large. Neither is it necessary to assume a time-scale sepa- 
ration between local and migration dynamics. We show that our results represent 
the typical transient of the metapopulation towards a quasi- steady state (if such a 
state exists). 

On the brink of extinction, that is, close to the bifurcation point where an in- 
finite metapopulation ceases to persist, we use a systematic expansion of an exact 
master equation in powers of A^~' (where is the number of patches) to find the 
most likely path to extinction, as well as the leading contribution to the time to ex- 
tinction. In this case (close to the bifurcation) the metapopulation dynamics sim- 
plifies considerably: it can be approximated by a sim ple, one-dimensional dynam- 
ics. T his fact is a consequence of a general principle (|Guckenheimer and HolmesL 



19831) stating that a dynamical system close to a bifurcation exhibits a 'slow 



mode': a particular linear combination of the dynamical variables is found to 
relax slowly, and the remaining degrees of freedom relax much more quickly and 
may be assumed to be in local equilibria. In other words, when the dynamical sys- 
tem is in a perturbed state, the slow mode evolves towards the equilibrium state 
on a longer time-scale than the fast variables do. This renders the dynamics effec- 
tively one-dimensional. We find the slow mode of the metapopulation dynamics 
and show how it depends on the properties of the local dynamics (given by the 
local growth rate and the local carrying capacity). The slow mode determines the 
stochastic dynamics of finite metapopulations as well as the deterministic dynam- 
ics of metapopulations consisting of infinitely many patches. In the latter case 
we find that the slow mode obeys Eq. ([T]). We derive how the parameters c and 
e depend on the parameters determining the life history of the local populations. 
We show, however, that the variable Q is in general not given by the fraction of 
occupied patches as envisaged by Levins. But it turns out that Q approaches the 
fraction of occupied patches in the limit of large carrying capacities and large local 
growth rates (this is the limit of time-scale separation mentioned above). 

In other words, we have derived Levins' model, Eq. ([T]), from a stochastic, 
individual-based model of a finite metapopulation. We find that Eq. ([T]) is still 
valid (close to the bifurcation) even when there is no time-scale separation be- 
tween the local and the migration dynamics. This is the consequence of the exis- 
tence of a slow mode. 
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Lande et alj (|l998h have suggested an elegant stochastic generalisation of Le- 



vins' model, in an attempt to compute how the local patch dynamics affects the 
properties of the quasi- steady state of the metapopulation (and the average time 
to extinction of this population). Their main idea is to connect the extinction and 
colonisation rates in Eq. ([T]) to local processes. The extinction rate is calculated as 
the inverse expected time to extinction of a single patch at the carrying capacity 
(which is determined self-consistently), and colonisation is defined as the rate of 
a single migrant arriving to an empty patch, seeding a population that grows to 
th e carrying ca pacity. This scheme is persuasive but not rigorous. The predictions 
of Lande et aP (1998) have, to our knowledge, never been tested by comparisons 
to results of simulations of stochastic, individual-based metapopulation models. 
Therefore it is important that our approach allows us to compute the rates of ex- 
tinction and colonisation from first principles. In the limit of large patch carrying 
capacities and close to the bifurcation we obtain expressions (exact in the limit 
we consider ) that a re very similar, but not identical, to the relations proposed by 



Lande et all (Il998h . 



Li summary, we characterise the stochastic dynamics of metapopulations on 
the brink of extinction. Using an expansion of the exact master equation describ- 
ing the stochastic dynamics with the number of patches as a large parameter, we 
identify a slow mode in the metapopulation dynamics, regardless of whether there 
is a time-scale separation between local and global dynamics or not. We show 
under which circumstances widely used metapopulation models provide accurate 
descriptions of metapopulation dynamics. 

The remainder of this paper is organised as follows. In section [2] we describe 
the individual-based stochastic metapopulation model investigated here. We sum- 
marise how our numerical experiments were performed and briefly describe how 
to represent the metapopulation dynamics in terms of a master equation, and how 
to expand this equation in powers of A^^'^ Our results are described in section 
m We first discuss the limit N ^ oo, demonstrate under which circumstances 
metapopulations persist in this limit, and analyse the metapopulation dynamics. 
Second, we turn to stochastic fluctuations of finite metapopulations, and sum- 
marise our results on fluctuations in the quasi-steady state, the most likely path 
to extinction as well as the average time to extinction. Section |4] contains our 
conclusions. Appendices A, B, and C summarise details of our calculations. 



5 



2. Methods 



In this section we define the stochastic, individual-based metapopulation model 
that is analysed in this paper. We describe how our numerical experiments are 
performed and derive a master equation, Eq. (fT4l) . that is the starting point for our 
mathematical analysis of metapopulation dynamics. In general it is not possible 
to solve this equation in closed form. We demonstrate how an approximate so- 
lution can be obtained by expanding the master equation using the number A'^ of 
patches as a large parameter. In the limit of ^ oo, the dynamics reduces to a 
deterministic model. The equilibri um prop erties of this mod el were analysed by 



Casagrandi and Gattol (119991 l2002h . and bv iNachmanI (l2000h 



2.1. Stochastic, individual-based metapopulation model 

The model consists of a population distributed amongst A'^ patches, as illus- 
trated in Fig. [U In each patch, the local population dynamics is a birth-death 
process with birth rates bf, and death rates J,. Here / denotes the population size 
in a given patch (and b^ = do = 0). To simplify the discussion we assume that the 
rates are the same for all patches, but more general cases can be treated within the 
approach described in this paper. In the following we illustrate our results for a 
particular choice of birth and death rates: 

bi = r i birth rate , (2) 

di = jii + {r - p)i^/K death rate . (3) 

The parameter r is the birth rate per individual, ft is the density-independent per 
capita mortality, and K determines the carrying capacity of a single patch. For 
simplicity, we take /i = 1 hereafter (this corresponds to measuring time in units of 
the expected life-time of individuals in the absence of density dependence). The 
parameters occurring in Eqs. (|2I3I) are listed in Table [B which summarises the 
notation used in this article. 

In addition to the local population dynamics, the number of individuals in each 
patch can change because some individuals emigrate from their patch to other 
patches, or be cause immigrants arrive from o ther patches. 



Following iHanski and GyllenbergI (|1993l) . the migration process is modelled 



as follows: individuals emigrate from a patch at rate m,, where / is the population 
size in the patch in question (mo = 0). If the individuals migrate independently 
and with constant rates, m, is proportional to z: 

m,- = m i emigration rate . (4) 
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More complex migration patterns can be incorporated. If for example individuals 
moved to avoid overcrowding, the emigration rate would be density dependent. 

The emigrants enter a common dispersal pool, containing the migrants from 
all patches that have not yet reached their target patch. Each migrant stays an ex- 
ponentially distributed time in the pool, with expected value 1 Irj, before reaching 
the target habitat, which is chosen with equal probability among all patches. This 
process is illustrated in Fig. [TJ Migration may fail if the individuals die before 
reaching the new habitat; this is modelled by the rate ^ of dying during dispersal. 
In practice, the probability of successful migration depends on the background 
mortality of the individuals, on the time the migrating individual spends in the 
dispersal pool, and on additional perils individuals may be exposed to during dis- 
persal (e.g. increased risk of predation due to lack of cover, etc.). In summary, 
if there are M migrants, (77 -I- individuals leave the dispersal pool per unit of 
time, and the rate of immigration to a given patch is / = rjM/N. 

2.2. Numerical experiments 

In the direct numerical simulations, the population evolves in the following 
manner. First, at any given time, the local rates, Eqs. dlHll), sum to a rate of the 
next locally generated event, A/t, for patch k. This rate is simply the sum of the 
local rates, Eqs. If patch k contains i individuals then we have: 

Kk = bi + di + mi . (5) 

The sum of these rates over all patches, A = X/t ^k, for = 1, . . . , A^, yields the 
rate for the next event occurring in the population. Thus the simulation proceeds 
by generating an exponentially distributed random number with expected value 
1/A. Second, at each time step, a patch is chosen with probability A^/A. Third, 
the type of event is chosen randomly: birth with probability b}l Ak, death with 
probability djlAk or emigration with mil A^. Fourth, the numerical experiments 
described below were performed in the limit 77 ^ 00 and for ^ = 0. This means 
that emigrating individuals are immediately assigned to a randomly chosen patch 
(possibly the one they come from). 

2.3. Master equation 
A natural and commonly adopted app roach to describe st ochastic population 



dynamics is to derive a master equation ( van Kampen , 198lb for the change in 
time of the probability p of observing ii individuals in the first patch, i2 individu- 
als in the second patch, . . . , and of observing M individuals in the migrant pool. 
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Recently this approach was adopted by iMeerson and Sasorovl (|201lh to describe 



the dynamics of local birth-death processes coupled by nearest-neighbour interac- 
tions (diffusion). 

In the following we pursue a different approach. Since the local population 
dynamics is assumed to be the same within all patches, it is sufficient to count the 
number of patches with a given number of individuals rather than keeping track 
of the number of individuals in each patch. Let nj denote the number of patches 
with j inhabitants at a given time. The state of the population is described by the 
variables no,ni, . . . , and M. In the master equation, the time derivative of the prob- 
ability p(«o, rii, M;t) of observing the system in a given state no,ni, . . . , M at 
time t, is the rate of arriving to a given state from other states, minus the rate of 
transitions to other states. Thus, we find the master equation for the probability 
p(nQ, Hi, . . . , M; t) by considering all possible transitions between the states of the 
population, contributing to the change of p(no, n^, . . . , M;t). 

For example, consider the effect of local births on the probability of finding 
the system in a given, 'focal' state no,ni, . . . , M at time t. A birth event in a patch 
with j individuals corresponds to the transition nj nj - I and nj+i nj+i + 1. 
Thus, dp/dt has the contribution 

-bjnjp(no,ni,...,M;t) (6) 

from births in the focal state. This contribution is negative since any such birth 
in the focal state moves the system away from this state. In order to obtain the 
positive contributions to dp/dt, we calculate the rate of arriving to the focal state 
by a birth in a different state: 

bjinj + l)p(. ..,nj + l, rij+i - 1, . . .) (7) 

A compact way of writing these contributions is in terms of the raising and low - 



ering operators E*, defined by their effects on a function g (Ivan Kampenlll98l|) 



E'^gi..,nj,...) = g(...,nj + \,...), (8) 
E']g(...,nj,...) = g{...,nj-l,...). 

In our example, the total contribution to dp/dt from birth events can thus be writ- 
ten as 

oo 

J](E'jE.,,-\)bjnjp. (9) 
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Adding up the contributions due to death, emigration, and immigration we find: 

. oo 

oo CO 

j=0 7=0 

OO 

+ ^ (EJEt^iE;, - l) 7jM% + (E; - 1) ^Mp . (10) 

7=0 

To simplify the discussion we assume that migration is instantaneous, this corre- 
sponds to taking the limit j] ^ oo and ^ = 0. In this limit, the immigration rate to 
a patch is given by: 



^=7^2"'^"^' ^^^^ 



7=0 

and the corresponding master equation takes the form: 

J CO oo 



7=0 7=0 

CO CO 



+ - Z Z E7_jEf E|ET^im,-/7,(ny - 6ij + 6i-ij)p - ^ nump . (12) 

1=1 7=0 !=1 

The last two terms on the right-hand side of Eq. (fT2l) describe instantaneous mi- 
gration where M = 0. The terms involving Kronecker 5-symbols (see Table [B on 
the right-hand side of Eq. (fT2)) arise from enforcing that, for the migration of an 
individual, emigration must precede immigration. These terms are of higher order 
in A^-^ 

The number of patches is given by = YjJ=q nj. The following discussion is 
simplified by making this constraint explicit in the master equation (fT2l) . This is 
achieved by considering no, the number of empty patches, to be be a function of 
nun2,...: 

oo 

no = N-J]nj. (13) 

7=1 
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Using Eq. (fT3l) we find the following master equation: 



dp(n, t) 



= Yji^-E.^i - \)bjnjp(n,t) + _^(E;Et_j - l)djnjp(n,t) 



OO OO 

;■=! j=l 

OO OO OO 

+ _^(E7 jE^'E^ - l)m,n;(l 5n)p{n, - J] m,n,p . (14) 

i=l /t=l i=\ 

Here the components /iy of the vector n = ini,n2, . . -Y denote the number of 
patches with j > 1 individuals, and E* = 1. Eq. (fT4l) describes the stochastic 
population dynamics of the metapopulation model considered here. 

In the next section we describe the approximate method of solving Eq. ([14]) 
adopted in the following. It corresponds to a systematic expansion of the master 
equation, using the number A^^ of patches in the population as a large parameter. 
This approach does not require that the number of individuals per patch is large, 
neither is it necessary to assume a time-scale separation between local population 
dynamics and migration. Our approach makes it possible to accurately describe 
the dynamics of metapopulations consisting of many small patches. In the limit 
of infinitely many patches, to lowest order in the expansion, a deterministic meta- 
population dynamics is obtained. Stochastic fluctuations in metapopulations with 
a large but finite number of patches are described by the leading order of the ex- 
pansion. 

2.4. Expansion of the master equation 

When N is large we expect that the probability p of observing a given dis- 
tribution of n changes only little when patches change in population size. It is 
important to emphasise that no assumption is made concerning the size of the 
changes to the population size of any single habitat. It is merely assumed that 
there are sufficiently many patches that they form a statistical ensemble. In this 
case it is perfectly possible to have big jumps in the population sizes of individual 
patches in the model without violating the assumption that the distribution of n is 
smooth, and to allow the local population dynamics to depend sensitively on the 
local patch population size when there are few individuals i n the patch — mod; 



elling for example the effect of abundance on mating success (|Saether et al.L 12004 



Melbourne and Hastingsll2008|) . 
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It is convenient to express the the components of the vector n in terms of the 
scaled frequencies fj = nj/N. The probability p of observing /i,/2, . . . is related 
to the probability p by 

p(f;t)=p(Nf;t). (15) 

In contrast to the distribution of rij, the distribution of fj is expected to become 
approximately independent of for a large number of patches. A change of a 
single count Uj leads only to a change of 1 /N in fj. Since the probability p is an 
approximately continuous function of fj in the limit of large values of A^, we seek 
an approximate solution of Eq. (fT4l) by expanding this equation in powers of N~^. 

In the limit of A^^ ^ oo we expect that stochastic fluctuations are negligible, so 
that the metapopulation dynamics becomes deterministic. It takes the form of a 
kinetic equation for /: 

^ = v(/). (16) 

at 

The right-hand side of this equation, v(/), depends upon the details of the model 
analysed. For the model described in Section 12. 1[ the form of v(/) is given in 
Eq. (fTSl) . In the limit of ^ oo, the question whether or not the metapopulation 
may persist is answered by finding the steady states /* of the system (fT6l) . given 
by v(/*) = 0. Persistence corresponds to the existence of a stable steady state with 
positive components f* > for some values of j. If, by contrast, the only stable 
steady state is /* = 0, then the metapopulation will definitely become extinct. The 
time to extinction is determined by Eq. (fT6l) . and by the initial conditions. 

For large (but finite) values of A^ the metapopulation fluctuates around the sta- 
ble steady states /* of Eq. (fT6l ). These fluctuations can be described by expanding 
the master equation (fT4l) to leading order in A^ ^ The stochastic fluctuations are 
expected to be small when A^ is large, but they are crucial to the metapopulation 
dynamics as mentioned in the Introduction. When the number of patches is finite, 
the metapopulation must eventually become extinct (the state f = is the only 
absorbing state of the master equation). When N is large, the time to extinction 
from a stable steady state of the deterministic dynamics is expected to be large. By 
analogy with standard large-deviation analysis, we expect that the time to extinc- 
tion increases exponentially with increasing A^, giving rise to a long-lived quasi- 
steady state. Below we analyse its properties, and estimate the time to extinction 
of the metapopulation. Our results are consistent with the above expectation, and 
we show that the time to extinction depends sensitively on the parameters of the 
model (r, K, m, and A^). 
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2.4.1. Deterministic dynamics in the limit ofN oo 

To simplify the notation, we drop the tilde in Eq. (fT5l) so that p(/, t) is the 
probability of observing, at time t, a fraction /i of patches with one individual, a 
fraction /2 of patches with two individuals, and so forth. In d eriving the expansio n 
of the master equation (fT4l) we use the approach described by Ivan Kampenl(ll98lh . 
Assuming that p is a smooth function of /, the action of the raising and lowering 
operators, Eq. ([8]), can be written as 



(17) 



Wj=tx^{±N-'df). 



The master equation (1141) is expanded as follows. We replace nj by Nfj in Eq. (1141) . 
insert Eq. (1T71) . and expand in powers of A^~'. Keeping only the lowest order in 
we arrive at an equation for p that corresponds to deterministic dynamics of 
the form (lT6l ). We find that the components of v(/) are given by: 

Vjif) = (bj-i + /)/;•-! + (dj^i + mj+Ofj+i (18) 
- (bj + I + dj + mj)fj for j > 1 , 

CO 

vi(/) = /(I - 2 fk) + (di + m2)/2 - (bi + I+di + mi)fi . 

k=l 

Here / = f^kfk is the rate of immigration into a given patch, correspond- 
ing to Eq. (ITTI) . Since / depends upon /, the deterministic dynamics (I16I18I) is 
nonlinear. We note that Eqs. (I16I18I) co rrespo nd to the metapo pulation model 
suggested by Casagrandi and Gatto^999 ) and Nachman ( 2000 ). Here we have 
derived it by a systematic expansion of the exact master equation in powers of 
where is the number of patches. Our derivation emphasises the fact that 
Eqs. (I16I18I) approximate the metapopulation dynamics by a deterministic equa- 
tion. Thi s approximation improves as becomes larger, and becomes exact as 
N ^ oo. lArrigonil (|2003[ ) has shown that this limit holds under quite general as- 



sumptions for the underlying stochastic model, namely that the time-evolution of 
the probability measure over the sta tes of the m odel conver ges to the determi nistic 
time evolution as A^^ ^ oo. Icasagrandi and Gattol (|l999 ). lNachman ( 2000|) . and 
others have studied how the stability of the steady states of the deterministic dy- 
namics depends upon the parameters of the model. However, Eqs. (I16I18I ) cannot 
be used to determine how stochastic population dynamics affects metapopulation 
persistence. In order to take the stochastic fluctuations into account, it is necessary 
to consider the next order in 1 /N. 
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2.4.2. Quasi-steady state distribution at finite but large values ofN 

Consider a stable steady state of the metapopulation in the limit of infinitely 
many patches, that is, a stable steady state /* of the deterministic dynamics (fT6l) . 

Metapopulations consisting of a finite number A'^ of patches exhibit random 
fluctuations caused by the random sequence of birth-, death-, and migration events. 
When the number of patches is large, these fluctuations are expected to be small, 
and the components fj of the vector / are expected to fluctuate closely around 
those of the stable steady state /*. If the fluctuations around /* are small, the 
metapopulation may persist for a very long time (but extinction of this finite met- 
apopulation is certain, as explained above.)- This situation is commonly referred 
to as a 'quasi-stable' steady state. We analyse its properties by expanding the 
master equation to leadi ng order in , using a standard method that is referred 



to as 'WKB analysis ' (IWilkinson et al.L |2007() . as the 'eikonal approximation' 



(iDykman et al. U 19941). or as the 'large -deviation principle' in the mathematical 
literature (Freid lin and Wentzelllll984|) . 

We briefly outline this method in the remainde r of this subsection. For a 



comprehensive description , the reader is referred to lAltland and SimonsI (|2010l) . 
Elgart and KamenevI (|2004|) . and lDykman et al.l (|1994l) . The method has been suc- 
cessfully used to describe fluctuations in finite biological popu lations, describing 



for example the spreading of epidemi cs (IDykman et al.l. 120081) or the r isk of ex 



tinction of biological populations (see lOvaskainen and MeersonI (|2010|) for a re- 
view). 

In the quasi- steady state we expect dp/d? a: and seek a solution of the master 
equation of the form 



P(/) * exp[-//5(/) -I- higher orders in ^] 



(19) 



The function S (/) is commonly referred to as the 'action' . It not only depends 
upon /, but also on the parameters of the problem (r, K, and m in our case). This 
dependence is not made explicit in Eq. (fT9l) . In the limit of a large number of 
patches, this function determines the form of the probability distribution p(/) in 
the quasi- steady state, and its sensitive dependence upon the parameters r, K, and 
m. 

As mentioned above, p(/) is expected to be strongly peaked at /* in the limit of 
large values of A^. In other words, the quasi-steady state distribution concentrates 
on the deterministic fixed point /* as the number of patches tends to infinity. It is 
therefore convenient to define the action function such that 5(/*) = 0. When the 
number of patches is large, we expect the distribution p(/) to be Gaussian in the 
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vicinity of the steady state. Correspondingly, S (/) is expected to be approximated 
by a quadratic function of 6f = f - f* . By contrast, non-Gaussian tails of this 
distribution, corresponding to large deviations of / from /*, reflect the particular 
properties of the extinction dynamics of the metapopulation. 

The form of the function S (/) is determined by inserting the ansatz (fT9l) into 
the master equation (fT4l) . making use of Eq. (fTVl ), and expanding in A'^"^ One finds 
a first-order partial differential equation for 5(/): 

= h(/,^1. (20) 



a/ 

In our case, for the master equation (fT4b , we obtain: 

CO oo 

mf,p) = ^(e"^^'-^^ - \)bjfj + YJ^^"'-'-''' - \)djfj 

CO CO 

+ Z - ^^"^'f'fj (21) 

,=1 j=l 

CO CO 

-Hj^(e^.-'-/'.-/"-l)m,/,(l-;^A). 

7=1 lc=l 

Here p = {puPi, ■ ■ 0^, Pj = dS /dfj for j > 1, and po = 0. 

The solution of Eqs. (12012 1|) is found by recognising th at Eq. (|20l) is a so-called 
'Hamilton- Jacobi equation' (Freidlin and Wentzelllll984t) . As a consequence, the 
action S (/) can be determined by solving the set of equations 

df dH dp dH 

These equations have the form of Hamilton's equations in classical mechanics 
with configuration- space variables /. The variables 

-I 

are therefore referred to as 'momenta'. An important property of Eq. ((22)) is that 
the function H{f(t),p(t)) remains constant along the 'trajectories' {f{t),p{t)), that 
is, along the solutions of Eq. (I22l) . 
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How is the form of S (/) determined by these solutions? The quasi-steady state 
solution of Eq. (fT4l) corresponds to solutions of Eq. (I22l) satisfying 



Pit) ^ 



as t ^ -oo , fit) f as ? ^ oo . 



and Hifit),pit)) = 0. (24) 

To every such solution corresponds an action obtained by integrating the 
momentum along the path if it), pit)): 



5(/) = Jdr/^. (25) 

Here p^ denotes the transpose of the vector p. 

The boundary conditions (|24)) are motivated as follows. Eq. (l20l) enforces 
H = 0. Further, observe that the stable steady state /* of the deterministic dy- 
namics (fT6l) corresponds to a steady state if*,0) of Eq. (|22l) . This can be seen by 
expanding the function Hif, p) to second order in p around p = 

Hif, p) = p'vif) + ^p^m)P +■■■■ (26) 

Here v(/) is given by Eq. (fTSl) . and the symmetric matrix D(/) has elements 
Dij = d^H/dpidpj. The elements are given in appendix lAl Eq. (|26l) shows that the 
dynamics for p = corresponds to the deterministic dynamics given in Eq. (fT6l) . 
The stability of the steady state if*,0) is determined by the eigenvalues of the 
matrix 

D 



The matrix A(/) has elements A,y = dvi/dfj evaluated at / = /*. It is the stability 
matrix of the stable steady state /* of the deterministic dynamics (fT6l) . Its ele- 
ments can be obtained from Eq. (IA.2I) in appendix |Al Similarly, D is evaluated at 
/*. Assuming that the steady state /* is stable, the eigenvalues of A must have 
negative real parts. In our case it turns out that the eigenvalues are in fact negative. 
We write > Ai > A2 > ■ ■ ■ . The eigenvalues of J occur in pairs A^, and -A^. The 
steady state if* ,0) of Eq. (|22l) is thus a saddle. In other words, stochastic fluc- 
tuations allow metapopulations consisting of a finite number of patches to escape 
from the steady state /* (that is stable in the limit A'^ ^ 00) to extinction (/ = 0). 
In general there are (infinitely) many s uch paths satisfying the boundary condi- 



tions (|24l) . iFreidlin and Wentzelll (|1984|) formulated a variational principle for the 
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most likely escape path: in the limit of large values of A^, the metapopulation goes 
extinct predominantly along this path. The quasi-steady state distribution reflects 
this property: configurations / along this pat h are assumed with higher pr oba- 
bility. According to the principle described by I preidlin and Wentzell ( 1984 ). the 



most likely escape path is the one with extremal action, Eq. (|25]) . 

The picture summarised above is schematically depicted in Fig.[2l for the case 
where the metapopulation persists in the limit of N ^ oo. Fig. 2a shows the f-p 
plane. As explained above, the point (f*,0) corresponds to the stable steady state 
of the deterministic dynamics, where the infinitely large metapopulation persists. 
The point (0, 0) corresponds to extinction, it is unstable. Solving Eq. (I22l) and 
inserting p = yields the deterministic dynamics (fT6l ). connecting these two fixed 
points. In the limit of N ^ oo, the metapopulation dynamics is constrained to 
the /-axis and must approach /*, as the arrow on the x-axis in Fig. indicates. 
In other words, in the situation depicted in Fig. ^ extinction never occurs in the 
limit ofN^ oo. 

In finite metapopulations, the situation is entirely diff"erent. Fig. [2^ illustrates 
that there is a path from /* to / = 0, reaching f = at the so-called 'fluctuational 
extinction point' (0,p*). Along this path, the momenta p{t) assume non-zero val- 
ues. These variables characterise the sensitivity of the finite metapopulation to 
stochastic fluctuations (p = in the deterministic limit ^ oo as mentioned 
above). The form of S(f) (Fig. [2b) is determined by evaluating Eq. (|25l) along 
this path, satisfying conditions (l24l) . In the limit o f large values o f A^, th e average 



time to extinction of the metapopulation scales as (IDykman et al.L Il994|) 



T,,,=Acxp[NS(f = 0)]. (28) 

The coefficient A may depend on A^, as well as r, ^,and m. In the argument of the 
exponent, only the leading A'^-dependence is explicit in Eq. (l28l) . The argument 
of the exponential determines the sensitive dependence of the time to extinction. 
In one-dimen sional problems with a single c omponent / (the case discussed in 



the review by lOvaskainen and Meersonl(|2010f) ). Eq. (|25l) shows that S (0) is given 



by the shaded area in Fig. [2^. In the case of our metapopulation model, by con- 
trast, the vector / has infinitely many components. It is therefore not possible, in 
general, to find the most likely path from /* to the fluctuational extinction point 
explicitly. In practice one may truncate the dynamics by only considering a fi- 
nite number of variables fj (up to a maximal value of jmax)- This is expected to 
be a good approximation when j'^ax is taken to be much larger than the carrying 
capacity K, since fj ^ for j » K. In our subsequent analysis of the problem 
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we make use of the fact that the dynamics si mplifies considerably in the vicinity 



of a bifurcation of th e deterministic dynamics (iGuckenheimer and Holmeslll983 



Dvkmanetal.Lll994l) . 



As mentioned above, near the steady state /*, the distribution (fT9l) is expected 
to be Gaussian. This corresponds to an action quadratic in 5/ = / - /*: 



(29) 



The matrix C, which is the covariance matrix of the distribution (fT9l) multiplied 
by N, can be obtained from the linearised dynamics, Eq. (I22l) . Using p = dS /df 
we have 

6p = C'^6f. (30) 

According to Eq. (|24l) . the dynamics must obey H{f(t),p(t)) = 0. It follows from 
Eq. (|26|) that the linearised dynamics satisfies this constraint provided 



AC + CA' + D = 0. 



(31) 



This equation determines the covariance matrix C in Eq. (|29|) in terms of the 
matrices A and D. 



3. Results and discussion 

In this section we summarise our results for the population dynamics of the 
metapopulation model described in Section lXTl using the expansion of the master 
equation outlined in Section [Z4l This section is divided into two parts. We first 
discuss the limit of infinitely many patches. Second, we analyse the most likely 
path to extinction in finite metapopulations. We also summarise our results for the 
average time to extinction of the metapopulation, and demonstrate that it depends 
sensitively on upon the parameters of the model (r, K, m, and A^^). 

3.1. Infinitely many patches 

It was shown in Section [Z4l that in the limit of infinitely many patches, the 
metapopulation dynamics is describ ed by the deterministic equat ion (fT6l), corre 



sponding to the model proposed by Casagr andi and Gattd (11999 ) and Nachman 



(12000). In this subsection we briefly summarise our results on the persistence and 
the relaxation behaviour of the metapopulation model introduced in section [211 
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3.1.1. Metapopulation persistence in the limit of infinitely many patches (N ^ oo) 
For sufficiently large emigration rates m, the deterministic dynamics (I16I18|) 
has two viable steady states. The state / = is unstable, and there is a second 
steady state /* given by 



fJ=foU 



+ r 



k=i 



dk + mk 



^"[r-pj ril*/r)nj+l)nj+l+z)' 

where Y{x) is the Gamma function, /* is the rate of immigration into a patch in the 
steady state, and z = pK/{r-p). Eq. (I32l) is most easily understood by recognising 
that the deterministic dynamics (I16I18I ) takes the form of a master equation for the 
probabilities f , that a patch is o ccupied by j indivi duals. Eq. ([321) is equivalent t o 



Eqs. (5a,b) in (lNachmanl . l200 0') and to Eqs. (5,6) in dCasag randi a nd GattoLl2002h 



In Eq. (I32l) . the factor is a normalisation factor (equal to the frequency of empty 
patches in the steady state) determined by the requirement that /g + Z^i /,* = 1- 
Note that the factors in the product in Eq. (|32l ) depend upon the rate /*. This gives 
rise to a self-consistency condition for the rate of immigration /* in the steady 
state: 



r = J]mjfJ. (33) 



The steady state (|32l) is stable provided the eigenvalues of A (this matrix is 
introduced in the previous section and its elements are given in appendix [A]) have 
negative real parts. This is the case when the steady-state immigration rate /* is 
larger than zero. The steady- state immigration rate is expected to decrease when 
the emigration rate decreases. Consider for instance decreasing m, defined in Eq. 
dH), while keeping all other parameters constant. As m approaches a critical value, 
mc, we observe that /* 0. At the same time all components of /* tend to zero, 
and /q 1. Below this critical point, that is for m < m^, the stable steady state 
ceases to exist (and f = turns stable). In the limit of infinitely many patches 
(A'^ ^ oo), the persistence condition 

m > m^ (34) 

ensures that a stable steady state exists, with f* > for some non-zero values of j. 
This persistence criterion for infinitely large metapopulations is equivalent to the 
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criteri on suggested by IChessonI (119841) and re-derived by ICasagrandi and Gatto 
(120021) . namely that the metapopulation persists provided that the expected num- 
ber of emigrants from a patch with initially one individual and into which immi- 
gration is excluded i s gr e ater th an unity. In a slightly different form this princi- 
ple is quoted by Ha nskil (Il998b . namely that the expected number of successful 
colonisations out of a given patch during its lifetime in an otherwise empty meta- 
population should be larger than unity. We emphasise that this criterion (or any of 
the equivalent criteria, such as (|34|) ) cannot be used to determine how stochastic 
fluctuations in finite metapopulations affect their persistence. 

The critical value nic is obtained by analysing the steady-state condition (|33l) 



for m close to irir. We write 



m = m(;(l -I- d) 



(35) 



and expand the steady-state immigration rate in powers of 6 (note that /* vanishes 
at 6 = 0): 



I* = Ii6 + h5' + .... 



(36) 



The constant /i is given in appendix |Bl Expanding Eq. (|33l) . we find to lowest 
order in 5 a condition for mc. 



k=2 



dk + nick 



(37) 



Fig. [3] shows how the solution mc of Eq. (l37l) depends upon the carrying capacity 
K for the model introduced in section 12.11 The critical migration rate is shown 
as a solid line in the left panel of Fig. [3l in the following referred to as the 'crit- 
ical line' . Above this line, a metapopulation consisting of an infinite number of 
patches persists. Expanding the condition (|37]) for large carrying capacities K we 
find that (see appendix 0: 



Mr 



1 



2nK 



exp 



-KW- 



r- 1 



(38) 



Fig. [3] also shows f* as a function of j for six stable steady states corresponding 
to different values of m and K. Results similar to those shown in th e six panels on 
the rig ht-hand side of Fig. [3] are giv en in Fig. 2 in (|Nachmanl . |2000|) . and in Fig. 2 
a-c in (ICasagrandi and Gattoll2002|) . 

One observes how the number of empty patches /q* tends towards unity as the 
critical line is approached. To leading order in 6 we have: 



f* = \-c,6. 



(39) 
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The constant ci is given in appendix |Bl Similarly, the expected number of indi- 
viduals per patch approaches zero: 



z 



^jfJ-C26 (40) 
j=i 

as 6 ^ 0. The constant C2 is given in appendix |Bl Fig. |4] shows the average num- 
ber of individuals per patch, and 1 - /q as a function of m. Shown are numerical 
solutions of the steady-state condition (|33l) . solid lines, of Eqs. (|39l ) and (|40l) , 
dashed lines, and of direct numerical simulations, symbols, of the metapopulation 
model described in section 12.21 The numerical experiments were performed for 
= 50, 100 and 1000 patches and averaged over an ensemble of different reali- 
sations. In principle, in a finite metapopulation no stable steady states exist with 
f* > for some j > 1. But the larger the number of patches, the larger the av- 
erage time to extinction is expected to be. = 1000, it turns out, is sufliciently 
large for the parameter values chosen that extinction did not occur during the sim- 
ulations. Consequently we observe good agreement between the direct numerical 
simulations and the numerical solution of the deterministic steady-state condition. 
When the number of patches is small, by contrast, this is no longer the case. In 
order to determine the quasi-steady state in such cases, the simulations must be 
run for a time large enough for the initial transient to die out, but shorter than the 
expected time to extinction. Stochastic realisations leading to extinction during 
the simulation time must be discarded. 

3.1.2. Relaxation dynamics in the limit of N oo and Levins' model 

We continue to discuss metapopulations with infinitely many patches. In this 
section we analyse how the metapopulation relaxes to the stable steady state /*. 
This question is important for two reasons. First, when the relaxation time is 
much smaller than the expected time to extinction (that is, when the number of 
patches is large enough), Eq. (fT6l) describes the relaxation dynamics well. Sec- 
ond, and more importantly, the deterministic dynamics exhibits a 'slow mode' in 
the vicinity of the critical line (in Fig. O: for small values of 6, it turns out, there 
is a particular linear combination of the variables fj that relaxes slowly, with a 
rate proportional to 6. This slow mode dominates the relaxation dynamics which 
is thus essentially one-dimensional f or small values of 6. This fa c t is w ell known 



in the theory of dyna mical systems (iGuckenheimer and HolmesL Il983b . see also 



Dykman et al.l (|1994l) . For our model, this fact has important implications. Es- 
sentially, the deterministic dynamics close to the critical line in Fig. |3]reduces to 
Levins' model, as we demonstrate in this section. 
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Fig. \5\ shows the relaxation of the average number of individuals and of the 
frequency of empty patches fo to their values in the stable steady state /*, for three 
different initial conditions. Shown are numerical solutions of Eqs. (I16I18I) . solid 
lines, as well as results of direct numerical simulations, symbols. The numerical 
experiments were performed for N = 50, 100 and 1000 patches, and averaged 
over 100 realisations. Even for 50 patches, and for the parameter values chosen 
in Fig. [51 the average time to extinction is large enough so that extinction did not 
occur in the numerical simulations. The agreement between the direct numerical 
simulations for = 1000 and the solution of Eqs. (I16I18I) is good. For N = 50 
and 100 we can observe deviations to the solution of Eqs. (|16I18I) . Here the effect 
of the finite number A'^ of patches becomes apparent. 

In general Eqs. ( I16I18I) must be solved numerically. We now show that the 
deterministic dynamics simplifies considerably close to the critical line in Fig.[3l 
when 6 is small. At 5 = we have that /* = 0, = 1 and consequently /* = 0. 
For small positive values of 6, the components of /* are of order 6, and we expand 



Eq. (fT6l ) in powers of 6 and / (|Dykman et al.L Il994|) 



= Z + "^^^ Z ^'fj + + 5) z ^Sm • (41) 



Here A|°^ are the elements of the matrix A evaluated at the critical line (6 = 0). 



'J 



They are obtained from (IA.2I ): 

Af = bi.,6i.ij + {du, + mS + l))6i^,j (42) 

- (bj + di + mj)6ij for i > I , 

Afj = (d2 + 2mc)6j2 - {bi +di+ m^)6ji + mJ . 

The slow mode exists because this matrix has one zero eigenvalue, All other 
eigenvalues are negative. At finite but small values of 5 the corresponding matrix 
(IA.2I ) evaluated at the steady state /* has one small eigenvalue, Ai ~ -5. To 
identify the slow mode we diagonalise A*^"^ given by (|42|) . Since this matrix is not 
symmetric, its left and right eigenvectors differ: 

A(°)/?, = ^f/?,, LjA<'« = Lji^°\ fora,;e= 1,2,.... (43) 

Here = (Li^,, L2a, ■ ■ •) denotes a left eigenvector with components Li^,, L2q, 

Ra denotes a right eigenvector with components Ria^Ria, ■ ■ ■■ We take the eigen- 
vectors to be bi-orthonormal, LlRa = Sap- As before, the eigenvalues are ordered 
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as > Af^ > Af^ > .... Multiplying Eq. (|4l]) from the left with and making 
use of Eq. (fT6l) yields the equations of motion of Qa = L^f (that is, the equations 
governing the time-evolution of Qa): 



dt 



(44) 



+ 



Aim = (that is, for 5 = 0) we have that Ai = = 0. For small values of 



6 we find \A\\ «; \Aa\ for a > \. While Aa, a > \, approach constants as 6 
Ai = L|A/?i oc -S tends to zero in this limit. This implies 



dt 



dt 



for a > I. 



0, 



(45) 



In the following, Qi is therefore referred t o as 'slow mode', w hereas Qa for a > 1 
are termed fast variables. As explained by lDykman et al.l (11994,) . the fast variables 
rapidly approach quasi- steady states 



Qa^- 



aa\5Q\ +aa2Ql 



i(0) 



for a > 1 



(46) 



that depend on the instantaneous value of the slow variable, Qi. Terms including 
fast variables (Qa for a > 1) are not kept on the right-hand side of Eq. (I46l) because 
they are of higher order in 6. This is a consequence of the fact that close to the 
steady state, Qi is of order 6. Eq. (|46l ) is correct to order 6^. The coeflicients a^i 
and Uai are given by 



aal = "^c ^ LaiA\^^Rij , 



(47) 



for or = 1,2, The elements Afj' and A|^^ are given in appendix lAl The equation 

of motion for the slow mode 2 1 is to third order in 6: 



dQi J 
— = au6Qi + anil + 5)Q\ 



dt 



(48) 



-I- terms of order 5^ involving Qp for jS > 1 
-I- terms of higher order in 5. 
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In the remainder of this subsection we neglect the 5^ -terms and higher-order terms 
in Eq. (|48] ). The 5^-terms are discussed in subsection 13.21 To second order in 6, 
the equation of motion for Qi is: 



(49) 



According to Eq. (I47I) . the coefficients an and an are given in terms of the compo- 
nents Lij of the left eigenvector Li and the components Rij of the right eigenvector 
/?i of A'^*'^ For the right eigenvector we find, from Eqs. (I42l) and (1431) . 



R 



1; 



k=2 



dk + ntck 



for J > 1 , 



(50) 



the first component being Rn. For the left eigenvector, we obtain the following 
recursion (with the boundary condition L^j = for j = -1): 



Lij+i 

This recursion is solved by: 



dj + Mci 
Lij = (Lij - Lij-i) 



-12 



-11 



1 -I- ^ I , and 



Ln 



(51) 



(52) 



Ln 



7-1 



bi by ^ nRu 



n=2 k=n+l 



for J > 2 . 



We show in appendix O that the elements of L| approach the following limiting 
form as ^ ^ oo: 



-11 



-12 



Z.11 



r + 1 



and L 



1; 



Ln r- 



\-r-J 
r- 1 



(53) 



We choose Ln such that Liy approaches unity for 7 » 1 as ^ ^ oo. This corre- 
sponds to the choice Ln = (r - l)/r, resulting in 



Li; = 1 



(54) 



in the limit of K ^ 00. In this limit, and for large values of r, we see that L| 
approaches the vector (1,1,...). The convention leading to Eq. (|54l) also fixes Rn 
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which must be chosen so that L^Ri = 1. Explicit expressions for an and a^, 
obtained from Eqs. (|47l) . (|50l) . and (|52l ) are given in appendix |Bl 

Eq. (|49l) has two steady states Qi = and 2* = -anS/an- Note that 2* is 
positive since an < 0. These two steady states correspond to the steady states 
/ = and /* of Eq. (fT6l) . Comparison with Eq. (I32l) shows that /?i is proportional 
to /* to lowest order in /* (that is, to lowest order in 6). In fact we have 

f* = RiQl + higher orders in 6 . (55) 

Fig. [6] illustrates how the deterministic dynamics relaxes to the stable steady 
state /*. Shown are the expected number of individuals per patch, YjJ jfj and 
the fraction of empty patches as functions of time, determined from numerical 
solutions of Eqs. (I16I18I) . Curves for three different initial conditions are shown 
(solid lines). The corresponding solutions of Eq. (|49l) are shown as dashed lines. 
Initially the variable Qi is not much slower than the Q^-variables for a > I, and 
the approximate one-dimensional dynamics (|49l) is not a good approximation. But 
the solution of Eqs. (I16I18I) rapidly relaxes to a form where Qi becomes a slow 
mode, and Eq. (|49l ) accurately describes the slow approach to the steady state. 

Comparing Eqs. (|49l ) and ([T]) shows that the slow mode Qi obeys Levins' 
equation. The results of this subsection allow us to compute e and c in terms of r, 
K, and 6. We find to order 6 

c - e = auS , (56) 
c = - ai2+ terms of order S . 

The contributions to c of order 6 can be computed explicitly using the formulae 
derived above. For the sake of brevity we do not specify these contributions here. 
In Sec. |3.2.3| we give an explicit formula valid in the limit of large K. 

We see that that Levins' model describes the dynamics of the metapopulation 
regardless of whether the patches are strongly coupled or not, provided the met- 
apopulation is sufliciently close to criticality (that is, close to the critical line in 
Fig. |3]). But we emphasise that in general the slow variable is not the fraction of 
occupied patches as envisaged by Levins. The variable Qi is given by L|/. Fig.|7] 
shows how the components Lij of Li depend upon j. When the local population 
dynamics is fast compared to the migration dynamics (for r = 1.5 and K = 50), 
the vector is approximately given by (1, 1, 1, . . .), see also Eq. (I54l) . In this case, 
Qi is approximately equal to the fraction Q of occupied patches. This is the limit 
of time-scale separation considered by Levins. We have thus derived the coef- 
ficients appearing in Eq. ([T]) from a stochastic, individual-based metapopulation 
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model defined in terms of tlie life history of its inhabitants (given by local birth 
and death rates, as well as the emigration rate). The asymptotic behaviours of an 
and ai2 for large values of K are given by (see appendixO: 



an ai2 



When the patches are strongly mixed by migration, then the interpretation of Qi 
is different. Fig. |7] shows that for r = 1.05 and K = 10, the components Lij are 
roughly proportional to j in the relevant range of j (where f* is not too small). In 
this case, therefore, the slow mode is interpreted as the average number of indi- 
viduals per patch, which in turn is proportional to the immigration rate, Eq. (fTTl) . 

In summary, in this section we have shown how the deterministic metapopula- 
tion dynamics simplifies in the vicinity of the bifurcation (critical line in Fig. [3]), 
that is, for small values of 6. In this limit, the deterministic dynamics is essen- 
tially one-dimensional, and of the same form as the deterministic e quation for the 
fraction of occupied patches originally suggested by Levins ( 19691) . Commonly it 



is argued that the form of Eq. ([T]) is appropriate when the local dynamics is much 
faster than migration. We have seen that this is not a necessary condition. In 
general, Eq. ([T]) describes the deterministic dynamics of metapopulations on the 
brink of extinction (for small values of 6). We find that the variable Q is indeed 
given by the fraction of occupied patches when the local-patch dynamics is fast. 
In general, however, Q has a different interpretation. Finally, we have been able to 
relate the rates c and e to the parameters of the individual-based, stochastic model 
(namely r, K, and 6 = {m- mc)/mc). 

Last but not least we emphasise that patch extinction in our model is entirely 
due to demographic fluctuations. It is possible to generalise our model so that 
the extinction rate a ccounts, in addition, for environrnental stochasticity, possibly 



including 'killing' (|Coolen-Schrijner and van Dooml |2006|) : transitions from an 



arbitrary state of a patch directly to extinction of that patch, local catastrophes in 
other words. 



3.2. Finite number of patches 

In the previous section we summarised our results on metapopulation dynam- 
ics in the limit of ^ oo. The subject of the present section is the dynamics for 
metapopulations consisting of a finite number A'^ of patches. In this case the fluc- 
tuations inherent in the birth-, death- and migration processes lead to fluctuations 
around the steady state /*. When the number of patches is large, these fluctua- 
tions are expected to be small. But they are essential: in a finite metapopulation 
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the only absorbing state is f = 0. In other words, the fluctuations turn the stable 
steady state /* of the deterministic dynamics into an unstable one. 

3.2.1. Fluctuations around the quasi-steady state 

For finite values of A'^, the population fluctuates around its quasi- steady state, 
as pointed out in section [T.4.2[ These fluctuations are Gaussian and of order A^~^ 
(cf. Eqs. (fT9l ) and (|29l)). Fig. [8^ shows how the the standard deviation of the 
fraction of empty patches, ctq, depends on the number of patches. A''. The standard 
deviation is calculated from the covariance matrix C, Eq. (|3TI) . as 

crl = if^) - {fof = -Yj ■ (58) 

Comparing the analytical approximation of ctq, calculated using Eqs. (I31I58I) . with 
simulations of the full stochastic model, we find good agreement when N > 100. 
Fig. [8h shows the variance for each cr^ = (ff) - (fj)^ as a function of j, for 
N = 100. 

3.2.2. Most likely path to extinction at finite but large values ofN 

When the number of patches is large, the state /* may be very long-lived 
(and is thus referred to as a quasi-steady state). The methods described in section 
12.41 allow us to systematically analyse the properties of this quasi-steady state in 
terms of Eq. (I22l) . While the deterministic dynamics in the limit of N ^ oo is 
given by Eq. (fT6l) . Eq. (|22l) describes the the stochastic fluctuations of the quasi- 
stable distribution at finite values of N prior to extinction. We note that by setting 
p = in the equations of motion [Eq. ((22l) ]. the deterministic dynamics [Eq. ([T6l) 1 
is obtained. In order to describe the fluctuations of a finite metapopulation around 
/* we need to find the most likely path from /* to a point / in the vicinity, namely 
the path from /* to / with extremal action [Eq. (|25])1. satisfying H{f{t),p{t)) = 0, 
as well as the boundary condition [Eq. (l30l)1. 

The tails of the quasi- steady state distribution near f = are obtained by 
computing the path from f* to f = with extremal action. This path is termed 
the most likely path to extinction, because in the limit of large values of the 
trajectories of the stochastic, individual-based dynamics to extinction are expected 
to fluctuate tightly around this path. The most likely path to extinction leaves the 
saddle point (/* , 0) along an unstable direction. Therefore this path cannot directly 
connect to the origin (0, 0). It turns out that the dynamics (|22|) exhibits a saddle 
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point dX f = and at non- vanishing momenta, which we call p* . The steady- 
state condition dpi At = gives rise to a recursion relation for the momentum 
components of p* of this saddle point: 

^Pl.rPl _ 1 = ^(1 _ qPI+pU) + _ QP'rPl^PU) . (59) 
bk bk 

The problem lies in finding the most likely path from (f*,0) to (0,p*). The man- 
ifold connecting these two points is infinite dimensional. It is straightforward 
to truncate the system of equations (f22l) at some large value of j, but finding 
the extremal path in the high-dimensional manifold (for instance by a numerical 
shooting method starting in the vicinity of /*) is very difficult. 

However, the problem simplifies considerably in vicinity of the critical line 
in the left panel of Fig. [3l When 6 is s mall, the dynamics (|22l) is approximately 



two-dimensional (IDykman et al.Lll994l) . since the linearisation (|27l) has two small 
eigenvalues when 6 <^ 1 . At 5 = we write 

j''^K = A^a^K, (60) 

J(0)<^' = -^(0)^' 

and corresponding equations for the left eigenvectors and H'J . Here J^"^ de- 
notes the matrix J (see Eq. (|27l) ) evaluated at 5 = 0, and Af^ are the eigenvalues 

of A'^K The left and right eigenvectors of J'^^' can be written in terms of those of 
A(0): 

K = ir]^ and K = i^A. (61) 



The left eigenvectors are given by 

£'J = (0, Rl), and £l = (lI O) . (62) 

The slow variables are obtained by projecting (f,py onto X| and X'l^. They 
are thus simply given by Qi and Pi, where 

Qa = Llf and P,=p''R,. (63) 



Following IDykman et al.l (119941) . the slow dynamics of Qi and Pi is found by 
expanding the Hamiltonian (|2TI) in powers of 6, f, and p. Starting from Eq. (I26l) . 
we use Eq. (|4T]) . as well as an expansion of D(/) in powers of /: noting that D 
vanishes at criticality, we write 

00 
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The elements of Bfjl = dDij/dfk are given in appendix |Al H{Qi,Pi) has the form 
(to third order in 6) 

HiQuPi) = PiianSQ, + anQ\) + bnQxP\ • (65) 

The coefficients 

bn = \y.LuB^'lUjR,,, (66) 

and an, an are given in appendixlB The equation of motion for Q\ and P\ is 

dQi dH dPi dH 

= and = -- — . (67) 

d? dPi dt dQi ^ 

The dynamics determined by Eqs. (I65I67I) is illustrated in Fig. |9l similar to Fig. [2l 
The three steady states of Eqs. (I65I67I) of interest for the questions addressed here 
are: 

2i = and Pi = , (68) 
21=0 and P = Pl = -anS/bu, (69) 
Qi = Ql = -and I an and P,=Q. (70) 

All three steady states are saddle points. The steady state (|70l ) corresponds to 
the saddle point (J* ,0) of the dynamics (|22l) . Note that an < 0. The steady 
state (|69l ) corresponds to the saddle point (0,p*) which lies at the end of the most 
likely path to extinction. In one-dimensional single-step birth-death processes, the 
corresponding point is commonly referred to as 'fluctuational extinction point' . To 
lowest order in 6 we find for the solution p* of Eq. (|59l ) 

p* = P*Li + higher orders in 5 . (71) 

Eq. dTT]) furnishes an interpretation of the left eigenvector L| for small values of 
6. This vector is proportional to the vector p*, and the components of this vector 
define the coordinates of the fluctuational extinction point. In one-dimensional 
birth-death processes, the corresponding value is given hy p* = - l o^Rg where Rn 



is the reproductive value. An ex ample is the so-called SIS-model (|Doering et al 



20051 : lAssaf and Meersonl 120101) . a stochastic model for the duration of the epi- 
demic state of an infectious disease. Here the reproductive value Rq is the expected 
number of infections caused, during its lifetime, by one infected individual intro- 
duced into a susceptible population. In infinite populations the epidemic persists 
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provided Rq > I. This criterion is precisely analogous to the persistence criteria 
discussed in Sec. I2.4.2I When Rq is only slightly larger than unity, Rq = \ + e, say, 

then p* « -e. 

Li our case, th e vector L| plays the role of a reproductive ygcfor (t pishei . 1930 ; 



SamuelsonL 119771) . Expanding the solution of the linearised deterministic dynam- 



ics (fT6l ), in terms of the eigenvectors of A*^°^ yields 

CO 

5f{t)= Yj^'"'R»Ll6m (72) 
~c^"Ri L|(5/(0) at large times . 



We see that the components Liy of the vector Li determine how much fluctuations 
of the number of patches with j individuals contribute to the relaxation towards 
the steady state in infinite metapopulations. In other words, the components of 
Li determine the susceptibility of patches with j individuals to stochastic fluctu- 
ations. We shall see below that this susceptibility determines the average time to 
extinction. We note that the correspondence (TTTI) implies that the solution of the 
recursion (|59l) for p* must be equivalent, to lowest order in 6, to the recursion (|5T1) 
for the components of Li, up to a factor. This is indeed the case. In analogy with 
the SIS-model discussed above, the coordinates of the vector L| parameterise the 
fluctuational extinction point. 

Finally we note that the components of Li have a simple interpretation in the 
limit of ^ ^ oo. Comparing Eqs. (|54l ) and (IC.19I ) we see that the component L^j is 
given by the probability of a single, isolated patch with j individuals to eventually 
reach its carrying capacity, K. This observation concludes our discussion of the 
nature of the fluctuational extinction point of Eq. (I67l) . 

The most likely escape path leads from the saddle (|70l) to the fluctuation ex- 
tinction point (|69l ). The path is parameterised by Qi and Pi. Solving H = for 
Pi we find from Eq. (|65l) (discarding the trivial solution Q\ = P\ = 0) 

Pi=-b\l{an5 + anQi). (73) 

This path corresponds to the straight line from (2], 0) to (0, P\) in Fig.|9l Eq. (1731) 
together with Eqs. (1551) and (TTTI ) imply that the /- and /^-spectra move rigidly 
towards extinction. In other words, our analysis shows that on the path to ex- 
tinction, both f{f) and p{t) retain their shape (initially given by the right and left 
eigenvectors R\ and L[ of A^"^) to lowest order in 5. 

In the limit of large carrying capacities (where the critical emigration rate is 
small) one might expect that patches become extinct independently of each other. 
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In this limit, the rate of extinction of single patches in the population is small, 
and the /-spectrum relaxes rapidly to its rigid shape once a given patch has gone 
extinct. It is important to emphasise that the patches are nevertheless coupled 
by migration which gives rise, in this limit, to a small rate of colonisation of 
empty patches. Thus the average time Text to extinction for the whole system 
is not determined by the largest time of extinction of the N single patches (we 
discuss the time T^^t in the following subsection). This has strong implications for 
conservation biology: it indicates the importance of protecting available (possibly 
empty) patches, even if they are small and prone to extinction. 

When the population is strongly mixed by migration, on the other hand, it is 
surprising that the /- and /^-spectra move rigidly to extinction. Migration upholds 
a balance that causes the shape of the /-spectrum to remain unchanged as the 
metapopulation comes closer to extinction. In other words, just the normalisation 
1 - /() changes. 

We have attempted to verify the predictions of this subsection by direct numer- 
ical simulations of the stochastic, individual-based model. The result is shown in 
Fig.[10l The simulations were performed as follows. For 10,000 independent re- 
alisations we followed the metapopulation to extinction. For each realisation we 
analysed the path to extinction by tracing the stochastic trajectory back in time, 
starting at the time of extinction. We traced the trajectories back for the average 
time it takes to reach, using this procedure, the number of individuals that cor- 
responds to the quasi-steady state. As a function of time, the average number of 
individuals per patch and the fraction of occupied patches was computed. Fig. 
[To! shows the probability of observing a given number of individuals in occupied 
patches, conditional on the fraction of occupied patches, 1 - /q. The probability 
is colour coded: white corresponds to high probability, black to low probability. 
The left panel of Fig. \T0\ is consistent with the prediction that the /-spectrum 
moves rigidly. This is not the case for the case shown in the right panel of Fig.fTOl 
Here 6 is too large, the two-dimensional approximation to the full Eq. (|22|) is not 
accurate. 

An alternative representation of the most likely path to extinction is depicted 
in Fig. \TT\ Employing the same simulations described above. Fig. [TT] shows the 
average number nj of patches with j individuals conditional on the total number of 
individuals, compared with results obtained by numerically integrating Eq. (|22|) . 
We observe that the agreement is good up to a point where nj becomes too small. 
When the total number of individuals is small (less than approximately 30 for the 
parameters in Fig.fTTI) discrete effects start to dominate, and the large- A/^ expansion 
of the master equation fails to accurately describe the last part of the trajectory 
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towards extinction. 



3.2.3. Time to extinction for finite but large values ofN 

The tail of the quasi-steady state distribution (fT9l) towards f = determines 
the average time to extinction, Eq. (|28]) : 

rext=Aexp(A^5). (74) 

The coefficient A may depend on A'^, as well as on r, K, and m. We have not been 
able to determine it. The action S = 5 (/ = 0) is a function of r, K, and m. Since 
NS appears in the argument of the exponential in Eq. (|74l ). the average time to 
extinction depends sensitively upon the number of patches, and on S . Close to 
the critical line in Fig. [3] we find the action by integrating Pi along the path given 
by Eq. (1731) . This yields: 

2 buai2 

which corresponds to the shaded area in Fig. |9l Note that ai2 < 0. The predic- 
tion for 5, Eq. (fTSl) . is compared to results of direct simulations in Fig. [121 The 
initial condition for the direct simulations was chosen as follows. First /* was 
calculated, this determines the initial state n = Nf* . Note however that the com- 
ponents nj of n must be integers. For small values of 6, all nj^o may round to zero, 
inconsistent with the constraint Yjj = N. In such cases we grouped counts nj 
corresponding to neighbouring values of j together before rounding. The action 
S was determined by plotting log Text as a function of A'^. For sufficiently large 
values of A'^ (such that NS » 1) one expects a straight line. We perform a linear 
regression by least squares to determine the slope S 16^. Fig. [12] shows log Text as 
a function of N6^. The approximations leading to Eq. ([751) require that A^ is large 
and 6 small. Eq. ([741) is expected to be a good approximation provided 

(A^5r'^2 « 5« 1 . (76) 

The data shown in Fig. [121 are consistent with this expectation. We see that the 
fitted values of S approach the analytical result, Eq. (TTSl ) as 6 is decreased (right 
panels in Fig. [121). This figure illustrates the sensitive dependence of the time to 
extinction upon r, K, 6, and A'^. 

In Fig. [131 we summarise our best estimates for S for different values of 
r and K (and for the smallest value of S for which we obtained reliable results). 
These estimates are compared to the analytical result, Eq. (TTSl ). We observe good 
agreement. Fig. [131 appears to indicate that S approaches the value 1/2 as ^ 
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increases. The approach appears to be the faster the large the value of r is. In 
order to demonstrate that this is in fact true, we have determined the asymptotic 
dependency of the coefficient bn upon K in the limit of large values of K. We find 
(see appendix 
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(77) 



This result, taken together with Eq. (l57l) determines the limiting value of S to 
be 1/2. Fig. [13] shows that already for r = 1.5 and ^ = 25 the value of S 
is very close to this limiting value. We note that while the expressions (|57l) and 
(TTTI) depend upon K, the action (fTSl) approaches a limit independent of K for large 
values of K. In this limit, the metapopulation dynamics can be understood in 
terms of a stochastic dynamics of the fr action of occupied p atches. 

Such an approach was suggested by Lande et al.l (|1998|) . In the following we 
briefly describe the similarit ies and differences b etween our asymptotic result and 
the approach suggested by Lande et al.l (|l998h . As already pointe d out in the 
introd uction, their analysis rests on four main assumptions: First, iLande et al 
(11998 ') consider the limit of fast local dynamics and slow migration (time-scale 
separation). I n our model, t his corresponds to the limit ^ ^ oo and to large values 
of r. Second, iLande et al. 1 ([1998) argue that the extinction rate in the stochastic 
model for the evolution of the number of occupied patches is given by the inverse 
time to extinction of a single patch, initially at carrying capacity. Third, it is 
assumed that the rate of successful colonisation is given by the product of the 
expected number of migrants and the probability that an em pty patch inva ded by 
one migrant grows to its (quasi-) steady state. Fourth, Lande et al. ( 1998h allow 
for the possibility that the number of individuals per patch in occupied patches 
changes as the metapopulation approaches extinction, and argue that this effect 
can be incorporated in terms of Q-dependent rates c{Q) and e{Q). 

Here we have shown, however, that in the limit of small values of 6 the dynam- 
ics of / on the most likely path to extinction is rigid: / does not change its shape, 
just its normalisation. This implies, in particular, that the average number of indi- 
viduals on occupied patches must remain unchanged as extinction is approached. 
This is clearly seen in Fig. [TOb . We note, however, that Fig. [TOb shows that for 
larger values of 6, f does change its shape (and the average number of individu- 
als in occupied patches decreases as extinction is approached). A first-principles 
theory for this effect is lacking. We refer to this point in the conclusions. 

Let us now con sider the parameterisations of the rates e and c suggested by 



Lande et al.l (|1998|) . Eqs. (3) and (4) in their paper. We have obtained analytical 
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results in the vicinity of the bifurcation (that is, for small va l ues of 6). In order to 
compare our results to the choices adopted by iLande et alj (Il998h we must take 
the limit K ^ oo. In this limit, it turns out, that the terms of order 6^ (and higher) 
in the second and third rows of Eq. (|48l) are negligible compared to the terms in 
the first row. In the limit of large values of K we therefore conclude: 



c = -(1 + 6)a 
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and 



e = c 



and . 



(78) 



The asymptotic expressions for an and au are given in Eq. (|57] ). In appendix 
O we have shown that these relations correspond precisely to the asymptotic K- 
dependence of the time Tk to extinction for a sin gle isolated patch a t carrying 
capacity. In the limit of large K (which is the limit iLande et al.l (|l998|) consider) 
we thus have: 

e=^. (79) 

This equation implies that the rate of extinction of patches in the coupled sys- 
tem does not depend upon 6 to leading order in 6. In other words, patches go 
extinct independen tly in the limit K ^ oo. Eq. (|79l ) closely resembles Eq. (4) in 



(LandeetaU 11998) 



Now consider the rate of successful colonisation. In the limit of ^ ^ oo, the 
probability of a single patch growing from one migrant to carrying capacity is 
U\K ~ (r - l)/r (see appendixjC])- Immigration is irrelevant in this context since 
the local patch dynamics is assumed to be much faster than migration. Using this 
expression for u\k and (|38l) for ntc we find the following expression for c 



c = nicil + 5)KuiK = rnKu\K 



(80) 



This result closely resembles Eq. (3) in (|Lande et al.L Il998h . We emphasise that 
our results are exact in the limit of small values of 6 and larg e values of A'^ and K. 

A further but minor diff"erence is that lLande et al.l (|1998b employ the diffusion 
approximation to evaluate their Eqs. (3) and (4). This approximation fails unless 
r is close to unity. 

We conclude this section with a discussion of Eq. (1751 ). Eqs. (I55I71I ) allow us 
to express the average time to extinction in the following form: 



log Text - ^{-p*^)r = y J]^-Pj)f; . 



N 
2 



(81) 



Here /* is the quasi-steady state distribution, and the vector -p*^ determines the 
susceptibility of patches with j individuals to stochastic fluctuations of the number 
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of individuals. As explained above, the components -p* are related to Fisher's 
reproductive vector characterising the susceptibility, for example, of age classes 
in matrix population-models. But note that in our case the components of -p* do 
not characterise the properties of classes of individuals, but of patches. In other 
words, Eq. (|8TI) can be understood by viewing the metapopulation as a population 
of patches, o r as a population of local populations, as originally envisaged by 



Levinsl ( 1969 ) 



3.2.4. The limit of large emigration rates 

Most results presented in this paper concern metapopulation dynamics close 
to the critical line shown in the left panel of Fig. |3l The reason is that the metapop- 
ulation dynamics simplifies substantially in this regime, as shown in the preceding 
sections. 

In this section we briefly discuss another limit of the metapopulation dynamics 
in the model introduced in Section 2, the limit of large emigration rates. This 
limit corresponds to the top part of the m-K plane shown in the left-most panel 
of Fig. [3l The limit is interesting for two reasons. First, in this limit too the 
metapopulation dynamics simplifies substantially. We show that the dynamics 
can be represented in terms of a one-dimensional process (just as Eq. ([T]) and the 
corresponding stochastic dynamics). But now the variable is the total number of 
individuals in the metapopulation and not Q\, or Q. Second, the results allow 
us to connect t he subject and the results of the present paper to recent results by 



HigginsI (120091) 



In the limit of infinite emigration rate (m oo) and no mortality during mi- 
gration = 0), the metapopulation behaves as a single large patch with modified 
birth and death rates that depend only on the total number of individuals in the 
metapopulation. This makes it possible to compute the average time to extinction 
exactly. This limit could describe, for example, biological populations in which 
larvae disperse randomly to all patches (e.g. pelagic marine species). 

We start from the individual-based, stochastic metapopulation model intro- 
duced in Section 2, in the limit 77 ^ oo and for ^ = 0. There are patches. It is 
convenient to characterise the state of the metapopulation by the number of indi- 
viduals in each patch /i , . . . , //v (rather than by the numbers /, denoting the fraction 
of patches with j individuals as in the preceding sections). Since migration is in- 
finitely faster than any other process, migration instantly brings the distribution 
of z'l, . . . , z'a? to multinomial form. Denoting the total number of individuals in the 
metapopulation by n, we see that after a birth (n - 1 ^ n) or death {n + \ ^ n). 
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the distribution instantaneously relaxes to the multinomial distribution: 



[ otherwise. 

We denote the birth and death rates for the whole metapopulation, conditional 
upon n, by B„ and D„, respectively. These rates can be calculated by averaging 
the local birth and death rates, Eqs. dH) and ©, conditional upon n. For the birth 
rate we find 

N N 

Bn = {J]bi,) = r{J^h) = m, (83) 

k=i k=i 

by virtue of Xf=i 4 = For the death rate we find: 

Dn= Yj (J]di^)Mnih,-..,iN)=^n + ^n^ + ^nil--\ . (84) 

i\,...,iN k=\ ^ ^ 

These rates (we take /i = 1 as before) define a one-dimensional, one- step birth- 
death process, that can be written in terms of a one-dimensional one-step master 
equation 

^ = (E- - 1 )Bn Pn + (E+ - 1 )Dn p„ , (85) 

at 

for the probability p„(0 of finding n individuals in the metapopulation at time t. 
The raising and lowering operators in Eq. ([85]) are defined in analogy with Eq. ([8]): 
E*g„ = ^,,+1. The average time to extinction for the one-dimensional Eq. (|85]) can 
be obtained by recursion (the result is given in appendix C). For the case of the 
multi-dimensional master equation (flOl) , by contrast, it was necessary to employ a 
large-A^ expansion before the dynamics could be reduced to the one-dimensional 
form ([U) on the brink of extinction. 



The results (I83H85I) allow us to investigate the question addressed by iHiggins 



(l2009h . namely how the risk of extinction depends upon the degree of 'fragmen 



tation' of a metapopulation, in the limit of infinite emigration rate. We define the 
total carrying capacity of the metapopulation ^tot = KN and ask how the aver- 
age time of extinction depends upon A'^, keeping ^tot constant. In other words, 
we fragment the metapopulation into N patches of equal size, such that the total 
carrying capacity remains the same. Figure [14] shows how the average time to 
extinction depends on the degree of fragmentation, for three different values of 
^tot (100, 500, and 1000), for r = 1.01 (other parameter values give qualitatively 
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similar results). The initial number of individuals is taken to be 5^tot (the aver- 
age time to extinction depends only weakly on the initial number of individuals, 
provided this number is larger than A'tot)- The analytical results (lines) are in good 
agreement with results, for large emigration rates, of direct numerical simulations 
of the model introduced in Section 2. Results are shown for two values of the 
emigration rate, m = 1 and 10. Each data point (symbols in Fig. [141) is estimated 
from 100 independent simulations. 

We find that the time to extinction decreases monotonically as the degree of 
fragmentation increases. For small degrees of fragmentation, the time to extinc- 
tion is approximately independent of the number of patches, but for larger levels 
of fragmentation, the time to extinction decreases markedly as fragmentation in- 
creases. Since the time to extinction is approximately exponentially distributed, 
this implies that the probability of extinction during a given time span increases 
monotonically as the metapopulation becomes more and more fragmented. 

This result is in qualitative agreement with the findings of Higgins (2009) for 
intermediate and la rge degrees of fragmentation. By contrast, when the number 
of patches is small, HigginsI (|2009() observes a decrease in the extinction risk with 
increasing fragmentation. A possible reason for this difference is that HigginsI 
fe009) considers local environmental fluctuations that would increase the local 
rate o f extinction compared to demographic fluctuations alone (Melb ourne and HastingsL 
l2008b . This could be important when the number of patches is low, but when the 
number of patches is large the fluctuations are averaged over and may become less 
important. 



4. Conclusions 

In this paper we have analysed metapopulation dynamics on the brink of ex- 
tinction in terms of a stochastic, individual-based metapopulation model consist- 
ing of a finite number A'^ of patches. In our model the distance from the bifurcation 
point where infinitely large metapopulations cease to persist is parameterised by 
the parameter 6. It measures the difference in emigration rate to its critical value 
at the bifurcation point. In the limit of large (but finite) values of and small 
values of 6 we have been able to quantitatively describe the stochastic metapopu- 
lation dynamics. We have shown that metapopulation dynamics, for small values 
of 6, is described by a one-dimensional equation of the form of Levins' model, 
Eq. ([U). We have derived explicit expressions for the parameters c and e in terms 
of the parameters describing the local population dynamics, and migration. We 
have shown that Levins' model is valid independently of whether or not there 
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is a time-scale separation between local and migration dynamics. The crucial 
condition is that the metapopulation is close to extinction. We note that in the 
absence of time-scale separation, the interpretation of the variable Q in Eq. ([T]) is 
no longer the fraction of occupied patches. More precisely, Q corresponds to the 
fraction of occupied patches if the growth rates and local carrying capacities are 
large enough; otherwise it has a different meaning, in particular, for small growth 
rates and carrying capacities it approximately corresponds to the average number 
of individuals per patch. 

The deter ministic limit of our model corre sponds to metapo pulation models 
discussed by I Casagrandi and Gatto ( 20021) and iNachman ( 2000 ). who have stud- 
ied the persistence of infinitely large metapopulations by analysing the stability of 
steady states. The corresponding persistence criteria (which we discuss in detail) 
do not allow to characterise the persistence of finite metapopulations. By con- 
trast, the stochastic dynamics derived and analysed in this paper makes it possible 
to characterise the most likely path to extinction and to estimate the average time 
to extinction of the metapopulation. We have discussed differences and similari- 
tie s between a particu lar asymptotic limit of our results, and the results obtained 
by lLande et al. (Il998 l). Fig.[l5] summarises different limiting cases of our results, 
and connections to earlier studies. 

While our approach results in a comprehensive description of metapopulation 
dynamics for the model we consider, many open questions remain. A technical 
point is that we have not yet been able to compute the pre-factor A in the expres- 
sion (|74l ) for the average time to extinction. This is a difficult problem requiring 
matching the solutions found in this paper with corresponding solutions valid for 
small values nj (the number of patches with j individuals). Such solutions are 
outside the scope of our large-A^ treatment. 

Many results derived in this paper are valid close to the critical line in Fig. [3l 
that is on the brink of extinction. The reason is that the metapopulation dynamics 
simplifies considerably in this regime, as we have shown in this paper. Further 
work is required in order to understand metapopulation dynamics for emigration 
rates much larger than the critical value. In this case, especially when the num- 
ber of migrants is large, we have observed deviations from Levins' model: the 
shape of the distribution of local patch population sizes is no longer rigid on 
the path to extinction. Biologically this reflects, for example, the rescue effect, 
whereby immigration from colonised patches extends the time to local extinction. 
A quantitative theory of this effect based on the life history of the individuals in 
the metapopulations is currently lacking. 

There is, however, a second limit of our individual-based, stochastic metapop- 
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ulation model that can be analysed in a simple fashion. This is the limit of very 
large emigration rates, m oo, where the metapopulation is well mixed, so that, 
in the absence of mortality during migration, it can be modelled as a sing l e larg e 
patch. In this limit it is possible to analyse the question posed by iHigginsI (|2009[) . 
namely how the risk of extinction of a metapopulation depends on the degree of its 
fragmentation. In the limit of large emigration rates, we have derived the expected 
time to extinction for our individual-based stochastic metapopulation model. We 
find that the time to extinction decreases as the degree of fragmentation of the 
population increases. There is excellent agreement between this prediction and 
results of direct numerical simulations of the model introduced in Section 2, for 
large but finite emigration rates. For hi ghly fragrnented populations, our con- 
clusions are consistent with the results of iHigginsI (l2009h . When the number of 
patches is small, by contrast, 'ffigginsj (t2009|) observes a different behaviour (dis- 
cussed in Sec. 13.2.41) . A possible reason for this difference is that IHigginsI (|2009|) 
considers local environmental fluctuations. These woul d increase the local rate of 
extinc tion compared to demographic fluctuations alone dMelboume and HastingsL 



20081) . This could be important when the number of patches is low. But when the 



number of patches is large the fluctuations are averaged over and are thus expected 
to be less important. To understand the differences and similarities between the 
models (and their biological implications) requires further work. 

The approach described in this paper allows us to address many other impor- 
tant questions in metapopulation dynamics. Consid er for example the dynamics of 
structured metapopulations. Core- satellite models dHanski and GyllenbergLll993|) 
of metapopulations have two types of patches, large and small. Large patches 
have more persistent populations, while populations in smaller patches are rela- 
tively ephemeral. Thus, the large patches are, on average, sources and the small 
patches are sinks for the metapopulation as a whole. Such models can be treated 
by suitable generalisations of the approach described in the present paper. 

Another important question is the role of environmental fluctuations. The 
within-patch population dynamics in our model takes demographic fluctuations 
into account, but does not explicitly incorporate the effect of local environmental 
fluctuations. Such fluctuations may be due t o, for example, yearly fluctuatio ns in 
growth rates, mortality, or carrying capacity ^ Melbourne and Hastings . 2008). The 
environmental fluctuations may be local, affecting patches independently. Alter- 
natively the fluctuations of the rates may be global, that is, the same for all patches 
(or at least highly correlated). 

What is the effect of local environmental fluctuations? Consider the set of 
patches with j individuals. One may expect that, when the number of patches is 



38 



large, only average rates enter into the deterministic equations of motion (fT6l) . But 
the environmental noise gives rise to fluctuations, and therefore leads to shorter 
persistence times of individual p atches. This in turn implies higher patch turnover 
jMelboume and HastingslboOSb . 

When the patches are few, the destabilising eff"ect of local env ironmental fluc - 
tuations on patches can significantly aff"ect the time to extinction (lHigginsLl2009h . 
The tendency of local environmental fluctuations to increase patch turnover could 
be represented in our model by introducing a separate process to the local popula- 
tion dynamics, whereby patches go extinct at a fixed rate, indep endently of the lo- 

cal population size, for example in terms of the 'killing' process (ICoolen-Schrijner and van Doom . 
l2006h . The theory outlined in this paper then makes it possible to determine the 
eff"ect of this process upon the time to extinction of the metapopulations. 

Global environmental fluctuations, affecting all patches simultaneously, can 
be incorporated in our model by introducing a time-dependent contribution 6fi{t) 
to the death rate in Eq. (3), the same for all patches. Taking this to be a random 
piecewise-co nstant function of ti me makes it possible to employ the approach 
described by ISchaper et al.l (|2012h . 

Last but not least, it is necessary to compare the predictions summarised here 
with those of models with an explicit spatial structure. We expect that the predic- 
tions summarised here should be the more accurate the wider the spatial scale of 
migration is. 
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Appendix A Formulae for the matrix elements occurring in the expansion 
of the master equation 

The matrix D(/) (see Eq. (|26l) ) has entries 
Dij = dij(ibi^i + 7)^1 + + m,+i)/;+i + (bi + I + di + mdfi) (A.l) 

+ - W+1 + ini+i)fi+i - {bi + /)/■) 

+ 6i.ij[ - (bi.i + I)fi.i - {di + mdf^ 

+ mi+ifi+i{fj_i - fj) + Mmjfj - nij^ifj+i) 

+ mifi{fj - fj^i) + fi-i{mj+ifj^i - nijfj) for /, ; > 1 , 

Dij = 6j2( - {d2 + m2)/2 - {bi + I)fi) + mj/sC/y-i - fj) + mji{fj - /y_i) 



+ Mmjfj-mj+ifj^i) + 



( oa \ 

V k=\ I 



{nij+ifj^i - nijfj) 



for J > 1 



Ai = Si2[ - {d2 + m2)/2 - {bi + I)fi) + m2f2{fi-i - fi) + mi/i(/- - /.i) 



+ Mrriifi - nii+ifi+i) + 



( oo ^ 

V k=\ J 



{mi+ifi+i - niifi) 



for / > 1 , 



Du =1 



+ {d2 + m2)/2 + {bi+I + di+ mi)/i 



V yt=l / 

/ DO 

+ 2 

. k=l 



i-Yjfk-fi 



(^2/2 -mi/i). 



The elements of A (see Eq. (|27l) ) are given by 

Aij = {bi-i + I)6i-ij + {di+i + mi+i)6i+ij - {bi + I + di + mi)6ij 
+ mj{fi-i - fi) 

CO 

Aij = {d2 + m2)6j2 - {bi + I + di + mi)6ji - I + mj(l - ^ A - /i) . 

k=i 



(A.2) 
for i > I , 



The elements of A'^^' in the expansion (|4TI) are obtained by setting 5 = in 
Eq. (IA.2I) . The elements of A^^^ occurring in the same expansion are found to 
be: 



i(i) 



j{6i+ij - 6ij) for i > 1 
i{5j2 - 5ji + 1) for ?■ = 1 



(A.3) 



40 



Finally, the coefficients A|^^, also appearing in the expansion (|4T)) . are given by: 

.(2) ^ f -K^ij - Si-i j) - j{5ik - 5i-]_k) for / > 1 
ijk I . + 1) _ + 1) for^=l • ^ • 

The coefficients (see Eq. (I64l) ) are given by: 

5^ = bu[5ki-x{5ij - 5ji.,) + 6i,{5ij - 5,7+i)] (A.5) 
+ (4 + m^)[6M+i{6ij - 5jM) + - 6ji.i)\ for z, j > 1, 

= ('^2 + m2){-6k25j2) - bi5k\5j2 + mk{5kj+\ - 5kj) for j > 1, 
^lil = (^<?2 + m2){-5k25i2) - bi6k\6i2 + mk{6ki+\ - Std for / > 1, 

= '^k + id2 + 3m2)42 + (^i + ^^i - mO^A-j. 

Appendix B Summary of results for coefficients used in section |3] 

In this appendix we give explicit expressions for the coefficients /i, ci, C2, an, 
an and bn appearing in Eqs. dlH), (gO]), (gg]), and dMl)- First, the coefficient 
/i in Eq. (I36l) is calculated by expanding Eq. (I33l) . Expanding to lowest order 
yields the condition (|37]) which gives mc but does not determine /i . In order to 
find /i it is necessary to expand Eq. (|33] ) to second order in 6. This requires 
expanding: 

°° ■'a + /* 

m,(i + 5) n j'l k = ~ ^2^' + ^i^^i^^' + ■ • • ' ^^-^^ 

j=2 k=2 " 



We find: 



j=2 



k 



di + mc "1 + '"c ^ -f?ii -fr^ + 
+ 



di+m,^~R^ ^b~i' 
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Inserting these expansions into the self-consistency condition (|33] ) gives: 

/, = ^. (B.4) 

^o(/i) 

We find 

Second, the coefficient ci determines the leading order of the ^-expansion of the 
fraction /(, of extinct patches in the steady state. To leading order in 6 we have 
f^{5) ~ 1 - ^o(^)- Thus the coefficient c\ in Eq. (|39l ) is given by 



7=2 



Third, the coefficient determining the lowest order of an expansion of the num- 
ber of individuals per patch, Eq. (I40l) . in powers of S is given by the lowest-order 
term of Eq. (IbH) . We find: 

C2 = — . (B.7) 

Fourth, in order to find the coefficients an and au appearing in Eq. (|49l ). we insert 
Eqs. (IA.4D and (IA.3I) into Eq. ([47]) for a = 1. Using Eqs. dSO]) and ^ yields 



aii=Liii?ii(ji-mcy— "^''•^ . ], (B.8) 



Finally, the coefficient bn in Eq. (|65] ) is determined by inserting Eq. (lA.ll) into 
Eq. (|66l ). which results in 

CO CO 

^11 = Yj^di + mJ)(Lu-Lu-ifRu-LnYj'^^^^^^'-^i'-^^^ii + ^^^nRn. (B.IO) 

(=2 1=2 
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Appendix C Asymptotics for large values of K 



In this appendix we briefly demonstrate how to obtain the asymptotics of the 
eigenvectors, of the coefficients an, a\2, bn, and of the critical emigration rate nic 
in the limit of ^ ^ oo (see Eqs. (|38] ) (|57] ). and (1771 ) in the main text). First, the 
results of this appendix demonstrate that the action (1751 ) tends to 1/2 as ^ — > oo. 
Second, the results obta i ned be low shed light on the connection of our results to 
the work of Lande etaP (Il998h . 

In the limit of large carrying capacities K, the asymptotics of the coefficients 
an, ai2, bu is given by the asymptotics of Rn which in turn is determined by the 
requirement that L|/?i = 1. The first step consists of deriving the asymptotic form 
of mc, Eq. (|38l) . in the limit of ^ ^ oo. The critical emigration rate nic is given by 
Eq. 07]): 



m. 



;=2 



k=2 



dk + nick 



(C.l) 



The product in this equation is estimated by exponentiation, taking the continuum 
limit, and the resulting integral is evaluated in the saddl e-point approximation , 
exact ly as describ ed by iDoering et al.l (|2005h . see also (iMehlig and Wilkinson , 
2007^) . Following Doering et al.l ( 2005 ) we introduce the variable x = i/K and 
define the functions b{x) and d{x) by 



bi = K b(x) and J, = K d{x) . 



(C.2) 



In the limit of ^ ^ oo, the term m^k in the denominator in (IC.ll) is negligible 
compared to dk- The product is then estimated using 



7-1 



1^ ^ j;j''d.vlogp(z)-i[logp(0)+logpCy)] ^ 



-K<5>(y) 



k=l 



(C.3) 



withy = j/K,p(x) = b{x)/d(x), and 

dx\ogp(x) = - - 
Jo ^ 



dxr 



+ (r - \)x 



(C.4) 



It follows that 



7=2 



k=l 



d, e-^-i-W 



(C.5) 
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The integral is evaluated in the saddle-point approximation. The saddle point is 
= 1. With d^O/dx^ = (r - l)/r at = I we have: 



y^fj^.^^ /^^±±eW-S)l, (C.6) 
^djl\dt \K{r-l)Kr^r 

resulting in Eq. (|38] ). 

The second step consists of determining the limiting form (|53] ) of the elements 
of Li , Eq. (|52l ). in the limit of large values of ^. We now show that the double sum 
in Eq. (|52l) grows exponentially, precisely cancelling the exponential decrease of 
mc as A" — > oo. The calculation is very c losely related to the e valuation of the time 
to extinction in a single-patch model by lDoering et al. (I2OO5I) . We use 



nRin Kr\K{r-\) Au^^^^ ^ ^ 

Performing the geometric sum and inserting the result into Eq. (l52l) we obtain 

L,j = \-r-L (C.8) 

The third step consists of estimating Rn which is determined by the require- 
ment 

We proceed as before and find: 

,„..^l^eW-5f)l. (CIO) 

These results enable us to determine the coeflicients an, ciu, ^12 from Eqs. (IB. 81) 
to (IB. 101 ). In Eq. (IB. 81) . the double sum is negligible compared to di in the limit of 
K ^ 00. This implies 

In Eq. (IB .91 ). the sum over j is negligible compare to unity. Evaluating the sum 
over k in the asymptotic limit, we find: 



an 



V 2n 
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Finally, Eq. (IB. 101) consists of three terms. The second term is negligible com- 
pared to the other two. This implies: 



V 2n 



(C.13) 



We see that the three coefficients have the same asymptotic dependence on K, 
apart from the minus sign in Eq. (IC.12I) . It turns out that this asymptotic K- 
dependence is exactly the inverse of the asymptotic dependence of the average 
time to extinction upon K for a single isolated patch with rates given by 
Eqs. (|2I3|) at carrying capacity K (no migration). Th e exact expression for the 



van Kampen. 


1981; 


Doering et al.. 


2005) 



^ n~ \ 1 I— I I 1 i-^ 'X J~ ^ 1 

^.=Ezn^EntEi-nl- 

" i=n+\ J 



i-\ 



n=\ k=\ 



(C.14) 



The correspo nding asympto t ic exp ression for / = K and large values of K is given 
in Eq. (19) of lOoering et aP (l2005h : 



Tk = 



In 



K(r- 1)3 



(C.15) 



Last but not least let us consider the probability u\k that a single i solated patch 
grows from one individual to its carrying capacity K. According to Ivan Kampen 
(Il98lb one has: 



U\K - 



di 



,=1 y=i '^n 



This expression becomes independent of K in the limit K 
same procedure as above, we find in this limit: 

r- 1 



UlK 



(C.16) 
oo. Following the 

(C.17) 



The c orresponding probability to grow from j individuals to K is (|van Kampen . 
198lh : 



UjK = 



(=1 k=\ 



bk 



K-l i , 



(C.18) 



1 + 



2.11 

(=1 k=l 
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This expression tends to 
as ^ ^ oo. 



UjK ~l-r 



(C.19) 
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Table 1 : Symbols used in this article 



Symbol Explanation 



Q Fraction of occupied patches, Eq. ([T]) 

c Colonisation rate, Eq. ([T]) 

e Extinction rate, Eq. ([T]) 

N Number of patches. Sec. 12.11 

r, Per-capita birth and death rates, Eqs. (I2I3I) 

K Patch carrying capacity, Eq. ^ 

m Per capita emigration rate, Eq. dH) 

bi, di, nii Birth-, death- and emigration rates for patch with / individuals, 
Eqs. (I213P1) 

/ Patch immigration rate. Sec. 12.11 

M Number of migrants in common migration pool. Sec. 12.11 

J] Rate at which migrants leave common migration pool. Sec. 12.11 

^ Rate of mortality during migration. Sec. 12.11 

Ayt Rate for next event generated in patch k. Sec. 12.21 

A Rate for next event occurring in population. Sec. 12.21 

n n = (ni,n2, . . .y where is number of patches with / individuals, 
Sec.O 

p{n, t) Probability of finding metapopulation in state n at time t, Eq. (flOl) 

Ej Raising and lowering operators, Eq. ([8]) 

5ij Kronecker delta: 6ij = 1 if / = j, and 6ij = Q\ii j. Sec. 12.31 

/ / = n/N: fi is fraction of patches with i individuals. Sec. 12.41 

p(f,t) p(fJ) = p(n,t), Eq. (fTSl ) (from Sec. 12.4.11 onwards, the tilde is 

dropped to simplify the notation) 

v(/) Right-hand side of deterministic equation for /, Eq. (fT6l) 

S{f) Action, Eq. (O 

p P = {puP2,...f: Pi = dSldfi,Scc.^A^ 

H{f,p) Hamiltonian function, Eq. (|2T]) 

A Jacobian matrix of the function v(/): A,y = dvi/dfj. Sec. 12.4.21 

D Matrix with elements Dij = d^H/dpidpj, Eq. ([261) 

J Jacobian matrix of the Hamiltonian dynamics, Eq. (ITTl) 

C Covariance matrix (multiplied by A'^): C,y = A^^ co\[fi, fj], Eq. (|3T]) 

0-5 Variance off/, for j > \: crj = Cjj/N, Sec lXTT] 

/* / at (quasi-)steady state, Eq. (I32l) 



Continued on next page. . . 
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Table 1 (continued): Symbols used in this article 
Symbol Explanation 

5f, dp Disturbance of / and p away from (J*,0), Eq. (I30l) 

/* Immigration rate at (quasi-)steady state, Eq. (I33l) 

p* p at the fluctuational extinction point, Eq. (|59l ) 

rtic Critical emigration rate, Eq. (|34l) 

6 6 = (m- mr)lmr. Eq. (|35l) 

A(°^ A evaluated at 5 = 0, Eq. (EB 

A^^^ Matrix with elements A^P = d^Vi/dfjdm, evaluated at 5 = 0, 

Eq. dH]) 

^5 ^SS = dhildfjdfudm, evaluated at 5 = 0, Eq. (EB 

Aa Eigenvalues of A, ordered such that A\ > A2> >■ • Eq. (|43l ) 

/if Eigenvalues of A(°), Eq. (|43]) 

Lq,, ^q, Left and right eigenvectors of A*^"^ with eigenvalues 1^"^ Eq. (|43]) 

J^°) J evaluated at 5 = 0, Eq. ([601) 

Ha, f^a Right eigenvectors of J^°^ with eigenvalues /if and -/if, 

Eq. (I6B 

Xff, Left eigenvectors of J*^°\ with eigenvalues /if and -/if, Eq. (I62l) 

Qa Slow (a = 1) and fast (or > 1) components of deterministic dy- 

namics for /, Eq. (|44l) 

Slow (or = 1) and fast (a > 1) momenta of the Hamiltonian dy- 
namics, Eq. (1631 ) 

55 5;]] = 5A7/5A,Sec.[l2j 

an, ai2, Coefficients in Eqs. (I48I65I ) 

Q\ Slow variable 2i at (quasi-) steady state, Eqs. (I55I70I ) 

P\ Slow variable P\ at fluctuational extinction point, Eqs. (16917 II) 

Text Expected time to extinction of the metapopulation, Eq. ([741) 

Tk Expected time to extinction of a single patch. Sec. 13.2.31 

UjK Probability that patch with j individuals reaches K individuals 

before becoming extinct, Eq. (|80l ) 
n Total number of individuals in metapopulation. Sec. 13.2.41 

Mn(i\, . . . , In) Multinomial distribution of z'l, . . . , i^, with 2f=i 4 = n, Eq. (I82l) 
5„, D„ Global birth and death rates in the limit m ^ 00, Eqs. (I83l[84l) 

Kiot Global carrying capacity in the limit m ^ 00, Sec. 13.2.41 
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Figures 




patch 1 patch 2 patch 3 patch 



Figure 1: Illustrates the stochastic, individual-based metapopulation model. The model describes 
local populations (also referred to as patches). The number of individuals in patch k is denoted 
by ik- Individuals are born and die with per-capita rates and di^. Furthermore, individuals 
may emigrate to a common dispersal pool (emigration rate m) where they stay an exponentially 
distributed time (with rate 77) before leaving the pool for one of the patches. Migration may fail 
if the individuals die before reaching the new patch, this possibihty is modelled by introducing a 
death rate ^ during dispersal. The instantaneous number of migrants in the pool is denoted by M. 
The innmigration rate from the dispersal pool into any given patch is given by / = r]M/N. 
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Figure 2: Illustrates the WKB method, (a) Schematic plot of the f-p plane. The significance of the 
variables / and p is explained in Section 1241 The deterministic dynamics, Eq. ( fTSI l. corresponds 
to motion along the /-axis. The form of the action Sif), Eq. dZSl l. is determined by the path 
from the quasi-steady state {f*,0) to the fluctuational extinction point {0,p*), corresponding to 
the most likely path to extinction. This is a standard situation, ofte n referred to in the literature. 
See for example Fig. 1 in Box 3 in ( Ovaskainen and Meerson , 2010l) . (b) Schematic. Shows 5'(/) 
as a function of / (solid line). Also shown (dashed line) is a quadratic approximation of 5'(/). 
This illustrates that the quasi-steady state distribution p(f) is Gaussian close to /*, but in general 
non-Gaussian in the tails. 
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Figure 3: Left: critical emigration rate Mc in the limit of infinitely many patches, — > oo, as 
a function of the carrying capacity K (solid line), computed from Eq. (l37b . The growth rate is 
r = 1.05. Above this critical line, metapopulations persist in the limit of — > oo. Right: six 
different stable steady states /*, obtained by solving Eqs. ( l32l[33] l numerically, for r = 1.05, and 
for the values of the emigration rate m and K indicated in the left panel. 
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Figure 4: (a) Average number of individuals per patch Yj°jLi jfj in the quasi-steady state as a 
function of the emigration rate m for three different values of the carrying capacity K. Shown 
are numerical solutions of Eqs. ( II6II8I 1, solid lines. Labels indicate the corresponding values 
of K. The growth rate is r - 1.05. Curves corresponding to the asymptotic expression ( |40] |. 
valid as m approaches the critical migration rate m^, are shown as dashed lines. Also shown are 
results of direct numerical simulations as described in section |Z2] for metapopulations consisting 
of N - 1000 (A), - 100 (□), and = 50 (O) patches, (b) Same but the fraction of occupied 
patches 1 -/^ as a function of m. The asymptotic behaviour as m — » nic is given by Eq. ( l39b (dashed 
lines). The results for the direct simulations are obtained as follows. For each set of parameters, 
30 stochastic runs up to a time smaller than the expected time to extinction are performed (if 
extinction occurs, the simulation is discarded). After the initial transient, 20 samples are taken 
from each simulation (separated by a time long enough so that the samples are uncorrected). If 
the simulations consistently become extinct during the initial transient we conclude that there is 
no quasi-steady state, and is set to unity. 
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Figure 5: (a) Average number of individuals per patch jfj as a function of time for r - 1 .05, 
K - 50, and m - 0.0167. Shown are numerical solutions of Eqs. ( 116118b for three different 
initial conditions: ten individuals per patch, 26 individuals per patch, and 50 individuals per patch 
as indicated by labels in the figure. Also shown are results of direct simulations as described in 
sectionlZ2l forA^ = 1000(a),A^ = 100(n)andA^ = 50(0), for each of the three initial conditions. 
The steady-state value is shown as a dashed line, (b) Same, but the frequency of empty patches 
/o as a function of time. Each point of the direct simulations corresponds to an average over 100 
stochastic realisations conditional on no extinction (no extinction actually occurred during these 
simulations). 
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Figure 6: Deterministic relaxation towards the stable steady state /*. (a) Average number of 
individuals per patch, jfj^ a function of time for r - 1 .05, K = 50, and 6 = 0.05. Shown as 
solid lines are numerical solutions of Eqs. ( I16I18I ) for three different initial conditions (indicated 
by labels in the figure): five individuals per patch, one individual per patch, and 10"^ individuals 
per patch on average. Also shown are solutions of Eqs. ( |46] | and ( |49] l (dashed lines), (b) Same, but 
the frequency of empty patches /o as a function of time. 
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Figure 7: Shows the components L\j of the left eigenvector L\ (see Eq. (143) )) as a function of j, 
for r - 1.5, K - 50, and for r - 1.05, K - 10 (solid lines). The dashed line corresponds to the 
limit A' — > oo, Eq. (|54] |. evaluated for r = 1.5. 
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Figure 8: Size of the fluctuations around the quasi-steady state of the distribution (T9[ . (a) Stan- 
dard deviation ctq of the fraction of empty patches as a function of the number of patches A^. 
Comparison of erg, calculated from C in Eq. ( |3TI ) and employing Eq. ( fSSl ) (solid line), to direct 
numerical simulations (O)- Parameters: r - 1.05, K = 50, m = 0.02. (b) Variances crj of the frac- 
tion of patches fj with j individuals as a function of j {j - 1,2,... ). Comparison of crj calculated 
using Eq. (ISTT i (solid line), to direct numerical simulations (O) for = 100. Parameters: r - 1.5, 
K^50,m^ 2.5 x 10 ^ 
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Figure 9: Illustrates dynamics in Qi-Pi plane, Eq. ( l67b . See Section [3.2.2l Similar to the schematic 
sketch shown in Fig.|2^. The path from (0, 0) to (Qf, 0) corresponds to the deterministic dynamics 
for Qi, Eq. ( |49] l. The path from the quasi-steady state to the fluctuational extinction point (Pj, 0) 
is a straight line. The action S , Eq. ( l75b . is given by the shaded area; S - Q\P\I2. 
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Figure 10: Results of direct numerical simulations of the stochastic, individual-based model, de- 
termining the most likely path to extinction. Shown is the probability distribution of the average 
number of individuals per occupied patch Yjj jfjl Tjjfj^ conditional on the fraction of occupied 
patches 1 - /q. Parameters: r = 1.5, K ^ 10, 6 ^ 0.103 (corresponding to m = 0.028), = 1000 
(a) and r = 1.05, K ^ 10, 6 ^ 0.61 (corresponding to m = 0.1342), and = 250 (b). The 
probability is colour-coded: high probability corresponds to white, low probability to black. How 
this plot was produced is described in Sec. 13.2.21 Also shown (only in panel a) is the result of the 
slow dynamics given by Eqs. (ISST l. (fTTT i. and (l73T l (solid green line). 
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Figure 11: Each panel shows the distribution of the number rij of patches with j individuals, 
conditional on that there are (n) individuals in the metapopulation. Thus the upper, left-most panel 
shows rij conditional on that there is one individual, and the lower right-most panel corresponds 
to 672 individuals (the expected number of individuals in the steady state). Results from direct 
simulations (averaged over 10"* stochastic realisations) (blue) are compared to results of numerical 
integrations of Eqs. ( 12 11221 ) (red). Parameters: r = 1.5, K^10,6^ 0.103 (m = 0.028), = 1000. 
Note that the first value of j in each panel is j=l. In other words, «(> is not shown. 
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Figure 12: (a) Average time to extinction Text from direct numerical simulations of the model 
described in section ITTl as a function of the number of patches A^, for r - \ .05, AT = 10. Symbols 
(X, O, A, V) correspond to the following values of 5: 0.0557,0.128, 0.182, 0.289, 0.611. Each 
data point corresponds to an average over 100 stochastic realisations. Also shown are numerical 
fits to the expected behaviour Eq. (l74b . (b) Action as a function of 5. Shown are the limiting 
result, Eq. dTSl l (dashed line), as well as the results of the fits shown in the left panel (O)- Error 
bars correspond to 95% confidence intervals. (C-d) Same as top panels, but for r = 1 .5 and K 
The symbols (x, Q, A, V, 0, <) correspond to 5 = 0.04, 0.064, 0.103, 0.143, 0.222, 0.301, 
and 0.458. (e-f) Same as top panels, but for r = 1.5 and K - 50. The symbols (x, O, A, V) 
correspond to 5 = 0.066, 0.173, 0.279, 0.493, and 0.706. 
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Figure 13: Action S divided by 6^ according to Eq. ( |75] |. in the limit of 5 —> 0, as a function of the 
carrying capacity K. Solid lines correspond to r = 1.05 and 1.5. Symbols (□, O, 0) correspond 
to the best estimates from Fig.[T2h. Fig.fTSb. and Fig.[T2fe. 
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Figure 14: Average time to extinction. Text, as a function of the degree of fragmentation Kiot/K. 
See Section [1211 Parameters: r = 1.01, m = 1, 10, Km = 100,500, 1000. Shown are results of 
Eqs. ( IC.14I83I84I I (solid lines) and results of direct numerical simulations as described in Sec. 13.2.41 
for m - 10 iO, □, 0) and m - I (V, A, <l). Error bars correspond to 95% confidence intervals. 
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Eq. (1), 

(Levins, 1969). 

This paper, section 3.1. 

n 

6^0 



Deterministic model for fj 
(Nachman, 2000; 
Casagrandi & Gatto, 1999). 
This paper, section 3.1. 



N ■ 



Eqs. (79,80). 

This paper, section 3.2.3. 

See also (Lande et al., 1998). 



K ->oo 
r s> 1 



Eqs. (65,67). 

This paper, section 3.2.2. 

5^0 



Stochastic model for fj. 



This paper, section 3.2. 



N —> oo 



Stochastic individual-based 
model, N patches. 
This paper, section 2.1. 



A' large but finite 



Instantaneous global dispersal. 

DO 

(Higgins, 2009). 

This paper, section 3.2.4. 



Figure 15: Schematic diagram showing relations between different limiting cases of the individual- 
based stochastic metapopulation model described in section ITTI The metapopulation consists of 
N patches coupled by migration. The variables fj denote the fraction of patches containing j 
individuals. The parameter m denotes the emigration rate, Eq. (|4|i, and 6 characterises how m 
differs from the critical rate Wc see Eqs. ( I35I37| | in Sec. 13. II 
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